Annoucements and such

install.packages(c("brms","tidybayes"))

Workspace setup

library(here)
library(tidyverse)
library(cowplot)
library(brms)
library(tidybayes)
library(patchwork)
# here's a bunch of stuff I include at the beginning of each qmd file. 
# it's mostly setting the color palette for figures, but it also controls the default code chunk option (showing you, i.e., echo) and setting the figure quality to ensure that it looks nice 
knitr::opts_chunk$set(fig.retina=3, echo=TRUE)
theme_set(theme_cowplot())
default_palettes <- list(
  c("#5e8485" , "#0f393a") ,
  c("#1c5253" , "#5e8485" , "#0f393a") , 
  # palette with 5 colours
 c( "#1c5253" , "#e07a5f", "#f2cc8f" , "#81b29a" , "#3d405b" ) ,
  # same palette interpolated to 8 colours
 c( "#1c5253" , "#e07a5f", "#f2cc8f" , "#81b29a" , "#3d405b" , "#a7a844" , "#69306d" ) 
  
)

options(ggplot2.discrete.fill = default_palettes, 
        ggplot2.discrete.colour = default_palettes)

Workflow

  1. State a clear question.

  2. Sketch your causal assumptions.

  3. Use the sketch to define a generative model.

  4. Use the generative model to build an estimator.

  5. Profit.


Model “recipes”

  1. Recognize a set of variables to work with. (Data and parameters.)
  2. Define each variable either in terms of the other variables OR in terms of a probability distribution.
  3. The combination of variables and their probability distributions defines a joint generative model that can be used to simulate hypothetical observations and analyze real ones.

Here’s an example:

\[\begin{align*} y_i &\sim \text{Normal}(\mu_i,\sigma) \\ \mu_i &= \beta x_i \\ \beta &\sim \text{Normal}(0,10) \\ \sigma &\sim \text{Exponential}(1) \\ x_i &\sim \text{Normal}(0,1) \\ \end{align*}\]


Model for globe-tossing

Here’s the model for last week’s globe-tossing experiment:

\[\begin{align*} W &\sim \text{Binomial}(N,p) \\ p &\sim \text{Uniform}(0,1) \\ \end{align*}\]

  • \(W\) is the observed count of water.
  • \(N\) is the total number of tosses.
  • \(p\) is the proportion of water on the globe.

The whole model can be read as:

The count \(W\) is distributed binomially with sample size \(N\) and probability \(p\). The prior for \(p\) is assumed to be uniform between 0 and 1.


Model for globe-tossing

Here’s the model for last week’s globe-tossing experiment:

\[\begin{align*} W &\sim \text{Binomial}(N,p) \\ p &\sim \text{Uniform}(0,1) \\ \end{align*}\]

Estimating the posterior using brms

Last week, we used grid approximation to estimate the posterior distribution. In the video you watched for today, McElreath moves on to something called QUADRATIC APPROXIMATION. It’s good to understand what that’s doing, but you and I are moving right along to MCMC. We won’t go into the details of what’s happening for a few weeks, but let’s start with the code to estimate the model.

Let’s say we tossed the globe 9 times and observed 6 waters:

m1 <-
  brm(data = list(w = 6),                             # <1>
      family = binomial(link = "identity"),           # <2>
      w | trials(9) ~ 0 + Intercept,                  # <3>
      prior(uniform(0, 1), class = b),                # <4>
      iter = 5000, warmup = 1000, seed = 3, chains=1, # <5>
      file = here("files/models/m21.1"))              # <6>
  1. Data can be a data frame or a list. Just make sure the variable names match what’s in your formula.
  2. How you assume your outcome variable is distributed.
  3. The formula for your outcome.
  4. Priors for every parameter in your model. Here, we only have one parameter, so we only need 1 prior. This happens to be a flat prior between 0 and 1.
  5. Some choices about how we want our model to run. We’ll go more into this later.
  6. These models can take a long time to run. You can automatically save the output to a file; when you do this, the next time you run this code, it won’t actually estimate the model, but will instead pull the output from your stated file. Be WARNED: if you change the data or the model code, it will NOT restimate your model until you delete the file.

Sampling from the posterior

Grid approximation gave us the calculated probability of each possible value of our parameter, \(p\). But our method of conducting bayes will no longer give us such a neat solution. Here’s how you get the posterior distribution for \(p\):

samples_from_post = as_draws_df(m1)
samples_from_post
## # A draws_df: 4000 iterations, 1 chains, and 3 variables
##    b_Intercept lprior lp__
## 1         0.46      0 -2.1
## 2         0.49      0 -1.9
## 3         0.56      0 -1.5
## 4         0.66      0 -1.3
## 5         0.71      0 -1.3
## 6         0.62      0 -1.3
## 7         0.68      0 -1.3
## 8         0.55      0 -1.5
## 9         0.63      0 -1.3
## 10        0.67      0 -1.3
## # ... with 3990 more draws
## # ... hidden reserved variables {'.chain', '.iteration', '.draw'}

samples_from_post %>%  
  ggplot(aes(x=b_Intercept)) +
  geom_density(fill = "grey", color = "white") +
  labs(x="Proportion water")


Posterior predictive distribution

ppd = posterior_predict(m1)
dim(ppd)
## [1] 4000    1
ppd
##         [,1]
##    [1,]    6
##    [2,]    4
##    [3,]    6
##    [4,]    5
##    [5,]    6
##    [6,]    6
##    [7,]    6
##    [8,]    4
##    [9,]    7
##   [10,]    7
##   [11,]    5
##   [12,]    7
##   [13,]    8
##   [14,]    7
##   [15,]    8
##   [16,]    9
##   [17,]    7
##   [18,]    9
##   [19,]    8
##   [20,]    9
##   [21,]    7
##   [22,]    9
##   [23,]    8
##   [24,]    7
##   [25,]    6
##   [26,]    6
##   [27,]    3
##   [28,]    4
##   [29,]    5
##   [30,]    3
##   [31,]    5
##   [32,]    8
##   [33,]    4
##   [34,]    7
##   [35,]    3
##   [36,]    2
##   [37,]    5
##   [38,]    5
##   [39,]    5
##   [40,]    5
##   [41,]    6
##   [42,]    8
##   [43,]    7
##   [44,]    6
##   [45,]    6
##   [46,]    6
##   [47,]    7
##   [48,]    6
##   [49,]    7
##   [50,]    9
##   [51,]    7
##   [52,]    3
##   [53,]    6
##   [54,]    7
##   [55,]    4
##   [56,]    7
##   [57,]    5
##   [58,]    2
##   [59,]    6
##   [60,]    6
##   [61,]    7
##   [62,]    6
##   [63,]    4
##   [64,]    6
##   [65,]    6
##   [66,]    6
##   [67,]    8
##   [68,]    7
##   [69,]    6
##   [70,]    4
##   [71,]    6
##   [72,]    5
##   [73,]    8
##   [74,]    5
##   [75,]    5
##   [76,]    3
##   [77,]    2
##   [78,]    6
##   [79,]    9
##   [80,]    5
##   [81,]    6
##   [82,]    7
##   [83,]    5
##   [84,]    5
##   [85,]    6
##   [86,]    7
##   [87,]    5
##   [88,]    5
##   [89,]    6
##   [90,]    6
##   [91,]    6
##   [92,]    3
##   [93,]    5
##   [94,]    5
##   [95,]    5
##   [96,]    6
##   [97,]    7
##   [98,]    7
##   [99,]    7
##  [100,]    8
##  [101,]    5
##  [102,]    7
##  [103,]    8
##  [104,]    6
##  [105,]    6
##  [106,]    8
##  [107,]    5
##  [108,]    6
##  [109,]    6
##  [110,]    3
##  [111,]    5
##  [112,]    4
##  [113,]    5
##  [114,]    8
##  [115,]    9
##  [116,]    3
##  [117,]    6
##  [118,]    5
##  [119,]    7
##  [120,]    7
##  [121,]    8
##  [122,]    3
##  [123,]    5
##  [124,]    7
##  [125,]    7
##  [126,]    7
##  [127,]    5
##  [128,]    6
##  [129,]    5
##  [130,]    6
##  [131,]    9
##  [132,]    8
##  [133,]    4
##  [134,]    5
##  [135,]    4
##  [136,]    8
##  [137,]    4
##  [138,]    7
##  [139,]    7
##  [140,]    5
##  [141,]    9
##  [142,]    5
##  [143,]    6
##  [144,]    5
##  [145,]    6
##  [146,]    8
##  [147,]    6
##  [148,]    6
##  [149,]    5
##  [150,]    7
##  [151,]    6
##  [152,]    9
##  [153,]    6
##  [154,]    7
##  [155,]    8
##  [156,]    4
##  [157,]    6
##  [158,]    7
##  [159,]    5
##  [160,]    3
##  [161,]    6
##  [162,]    3
##  [163,]    5
##  [164,]    8
##  [165,]    6
##  [166,]    6
##  [167,]    3
##  [168,]    6
##  [169,]    7
##  [170,]    9
##  [171,]    5
##  [172,]    5
##  [173,]    5
##  [174,]    8
##  [175,]    7
##  [176,]    3
##  [177,]    9
##  [178,]    9
##  [179,]    6
##  [180,]    7
##  [181,]    5
##  [182,]    7
##  [183,]    9
##  [184,]    7
##  [185,]    6
##  [186,]    5
##  [187,]    8
##  [188,]    6
##  [189,]    6
##  [190,]    5
##  [191,]    8
##  [192,]    8
##  [193,]    6
##  [194,]    7
##  [195,]    4
##  [196,]    5
##  [197,]    6
##  [198,]    7
##  [199,]    8
##  [200,]    8
##  [201,]    9
##  [202,]    4
##  [203,]    8
##  [204,]    7
##  [205,]    3
##  [206,]    6
##  [207,]    4
##  [208,]    6
##  [209,]    7
##  [210,]    7
##  [211,]    7
##  [212,]    4
##  [213,]    9
##  [214,]    5
##  [215,]    8
##  [216,]    8
##  [217,]    9
##  [218,]    6
##  [219,]    7
##  [220,]    5
##  [221,]    6
##  [222,]    4
##  [223,]    3
##  [224,]    6
##  [225,]    4
##  [226,]    3
##  [227,]    4
##  [228,]    9
##  [229,]    6
##  [230,]    6
##  [231,]    8
##  [232,]    7
##  [233,]    5
##  [234,]    5
##  [235,]    7
##  [236,]    6
##  [237,]    9
##  [238,]    7
##  [239,]    5
##  [240,]    3
##  [241,]    3
##  [242,]    6
##  [243,]    7
##  [244,]    8
##  [245,]    6
##  [246,]    6
##  [247,]    1
##  [248,]    8
##  [249,]    4
##  [250,]    5
##  [251,]    6
##  [252,]    8
##  [253,]    7
##  [254,]    9
##  [255,]    8
##  [256,]    7
##  [257,]    4
##  [258,]    7
##  [259,]    7
##  [260,]    4
##  [261,]    7
##  [262,]    8
##  [263,]    5
##  [264,]    8
##  [265,]    4
##  [266,]    8
##  [267,]    5
##  [268,]    7
##  [269,]    6
##  [270,]    8
##  [271,]    6
##  [272,]    6
##  [273,]    9
##  [274,]    5
##  [275,]    5
##  [276,]    4
##  [277,]    9
##  [278,]    8
##  [279,]    7
##  [280,]    7
##  [281,]    6
##  [282,]    3
##  [283,]    5
##  [284,]    1
##  [285,]    2
##  [286,]    1
##  [287,]    1
##  [288,]    9
##  [289,]    6
##  [290,]    8
##  [291,]    7
##  [292,]    9
##  [293,]    9
##  [294,]    7
##  [295,]    4
##  [296,]    6
##  [297,]    8
##  [298,]    6
##  [299,]    1
##  [300,]    2
##  [301,]    5
##  [302,]    7
##  [303,]    6
##  [304,]    7
##  [305,]    7
##  [306,]    6
##  [307,]    8
##  [308,]    7
##  [309,]    6
##  [310,]    5
##  [311,]    4
##  [312,]    7
##  [313,]    7
##  [314,]    4
##  [315,]    6
##  [316,]    8
##  [317,]    5
##  [318,]    7
##  [319,]    5
##  [320,]    4
##  [321,]    4
##  [322,]    4
##  [323,]    4
##  [324,]    3
##  [325,]    5
##  [326,]    7
##  [327,]    5
##  [328,]    5
##  [329,]    4
##  [330,]    8
##  [331,]    4
##  [332,]    7
##  [333,]    8
##  [334,]    5
##  [335,]    8
##  [336,]    8
##  [337,]    8
##  [338,]    5
##  [339,]    2
##  [340,]    3
##  [341,]    5
##  [342,]    6
##  [343,]    7
##  [344,]    7
##  [345,]    7
##  [346,]    4
##  [347,]    8
##  [348,]    9
##  [349,]    8
##  [350,]    7
##  [351,]    8
##  [352,]    6
##  [353,]    4
##  [354,]    4
##  [355,]    8
##  [356,]    7
##  [357,]    3
##  [358,]    1
##  [359,]    3
##  [360,]    6
##  [361,]    4
##  [362,]    4
##  [363,]    7
##  [364,]    4
##  [365,]    6
##  [366,]    4
##  [367,]    7
##  [368,]    4
##  [369,]    9
##  [370,]    8
##  [371,]    6
##  [372,]    3
##  [373,]    3
##  [374,]    3
##  [375,]    2
##  [376,]    3
##  [377,]    4
##  [378,]    8
##  [379,]    2
##  [380,]    6
##  [381,]    9
##  [382,]    6
##  [383,]    9
##  [384,]    7
##  [385,]    7
##  [386,]    7
##  [387,]    6
##  [388,]    5
##  [389,]    3
##  [390,]    7
##  [391,]    7
##  [392,]    7
##  [393,]    7
##  [394,]    4
##  [395,]    8
##  [396,]    9
##  [397,]    5
##  [398,]    5
##  [399,]    8
##  [400,]    7
##  [401,]    6
##  [402,]    6
##  [403,]    7
##  [404,]    8
##  [405,]    4
##  [406,]    6
##  [407,]    6
##  [408,]    7
##  [409,]    9
##  [410,]    8
##  [411,]    8
##  [412,]    7
##  [413,]    6
##  [414,]    4
##  [415,]    7
##  [416,]    6
##  [417,]    4
##  [418,]    4
##  [419,]    7
##  [420,]    4
##  [421,]    5
##  [422,]    7
##  [423,]    2
##  [424,]    6
##  [425,]    8
##  [426,]    6
##  [427,]    3
##  [428,]    8
##  [429,]    6
##  [430,]    5
##  [431,]    6
##  [432,]    3
##  [433,]    4
##  [434,]    7
##  [435,]    7
##  [436,]    8
##  [437,]    7
##  [438,]    3
##  [439,]    3
##  [440,]    8
##  [441,]    6
##  [442,]    8
##  [443,]    7
##  [444,]    7
##  [445,]    7
##  [446,]    4
##  [447,]    5
##  [448,]    5
##  [449,]    4
##  [450,]    7
##  [451,]    3
##  [452,]    5
##  [453,]    6
##  [454,]    2
##  [455,]    6
##  [456,]    3
##  [457,]    8
##  [458,]    7
##  [459,]    9
##  [460,]    6
##  [461,]    8
##  [462,]    7
##  [463,]    8
##  [464,]    5
##  [465,]    8
##  [466,]    7
##  [467,]    5
##  [468,]    4
##  [469,]    6
##  [470,]    5
##  [471,]    5
##  [472,]    8
##  [473,]    6
##  [474,]    5
##  [475,]    6
##  [476,]    7
##  [477,]    4
##  [478,]    5
##  [479,]    7
##  [480,]    6
##  [481,]    7
##  [482,]    5
##  [483,]    6
##  [484,]    6
##  [485,]    5
##  [486,]    6
##  [487,]    5
##  [488,]    4
##  [489,]    6
##  [490,]    6
##  [491,]    8
##  [492,]    4
##  [493,]    2
##  [494,]    3
##  [495,]    5
##  [496,]    7
##  [497,]    6
##  [498,]    5
##  [499,]    3
##  [500,]    3
##  [501,]    4
##  [502,]    5
##  [503,]    7
##  [504,]    5
##  [505,]    4
##  [506,]    6
##  [507,]    8
##  [508,]    6
##  [509,]    6
##  [510,]    7
##  [511,]    8
##  [512,]    6
##  [513,]    6
##  [514,]    6
##  [515,]    6
##  [516,]    6
##  [517,]    6
##  [518,]    7
##  [519,]    4
##  [520,]    4
##  [521,]    2
##  [522,]    1
##  [523,]    3
##  [524,]    5
##  [525,]    4
##  [526,]    8
##  [527,]    8
##  [528,]    5
##  [529,]    4
##  [530,]    7
##  [531,]    7
##  [532,]    9
##  [533,]    8
##  [534,]    5
##  [535,]    6
##  [536,]    5
##  [537,]    6
##  [538,]    7
##  [539,]    6
##  [540,]    7
##  [541,]    3
##  [542,]    5
##  [543,]    6
##  [544,]    7
##  [545,]    7
##  [546,]    8
##  [547,]    5
##  [548,]    6
##  [549,]    6
##  [550,]    4
##  [551,]    8
##  [552,]    7
##  [553,]    3
##  [554,]    6
##  [555,]    6
##  [556,]    4
##  [557,]    3
##  [558,]    3
##  [559,]    3
##  [560,]    4
##  [561,]    4
##  [562,]    6
##  [563,]    6
##  [564,]    5
##  [565,]    6
##  [566,]    4
##  [567,]    5
##  [568,]    7
##  [569,]    8
##  [570,]    7
##  [571,]    7
##  [572,]    2
##  [573,]    4
##  [574,]    4
##  [575,]    5
##  [576,]    8
##  [577,]    6
##  [578,]    7
##  [579,]    5
##  [580,]    6
##  [581,]    5
##  [582,]    5
##  [583,]    6
##  [584,]    6
##  [585,]    6
##  [586,]    5
##  [587,]    2
##  [588,]    5
##  [589,]    5
##  [590,]    2
##  [591,]    5
##  [592,]    5
##  [593,]    7
##  [594,]    6
##  [595,]    4
##  [596,]    6
##  [597,]    6
##  [598,]    4
##  [599,]    5
##  [600,]    8
##  [601,]    7
##  [602,]    4
##  [603,]    5
##  [604,]    7
##  [605,]    6
##  [606,]    8
##  [607,]    8
##  [608,]    2
##  [609,]    8
##  [610,]    6
##  [611,]    2
##  [612,]    9
##  [613,]    5
##  [614,]    5
##  [615,]    5
##  [616,]    3
##  [617,]    4
##  [618,]    8
##  [619,]    8
##  [620,]    7
##  [621,]    5
##  [622,]    7
##  [623,]    6
##  [624,]    7
##  [625,]    3
##  [626,]    4
##  [627,]    6
##  [628,]    9
##  [629,]    4
##  [630,]    4
##  [631,]    3
##  [632,]    7
##  [633,]    8
##  [634,]    5
##  [635,]    5
##  [636,]    5
##  [637,]    6
##  [638,]    0
##  [639,]    0
##  [640,]    3
##  [641,]    5
##  [642,]    5
##  [643,]    6
##  [644,]    7
##  [645,]    5
##  [646,]    5
##  [647,]    3
##  [648,]    4
##  [649,]    7
##  [650,]    5
##  [651,]    8
##  [652,]    6
##  [653,]    7
##  [654,]    2
##  [655,]    4
##  [656,]    7
##  [657,]    9
##  [658,]    8
##  [659,]    6
##  [660,]    8
##  [661,]    7
##  [662,]    5
##  [663,]    4
##  [664,]    5
##  [665,]    8
##  [666,]    7
##  [667,]    9
##  [668,]    8
##  [669,]    7
##  [670,]    5
##  [671,]    7
##  [672,]    7
##  [673,]    4
##  [674,]    4
##  [675,]    6
##  [676,]    8
##  [677,]    7
##  [678,]    5
##  [679,]    7
##  [680,]    3
##  [681,]    5
##  [682,]    6
##  [683,]    5
##  [684,]    7
##  [685,]    4
##  [686,]    7
##  [687,]    6
##  [688,]    5
##  [689,]    8
##  [690,]    5
##  [691,]    2
##  [692,]    2
##  [693,]    1
##  [694,]    5
##  [695,]    6
##  [696,]    8
##  [697,]    6
##  [698,]    2
##  [699,]    3
##  [700,]    4
##  [701,]    6
##  [702,]    5
##  [703,]    6
##  [704,]    3
##  [705,]    6
##  [706,]    3
##  [707,]    5
##  [708,]    6
##  [709,]    8
##  [710,]    8
##  [711,]    4
##  [712,]    7
##  [713,]    7
##  [714,]    5
##  [715,]    6
##  [716,]    7
##  [717,]    7
##  [718,]    4
##  [719,]    7
##  [720,]    3
##  [721,]    4
##  [722,]    2
##  [723,]    5
##  [724,]    5
##  [725,]    6
##  [726,]    5
##  [727,]    4
##  [728,]    4
##  [729,]    4
##  [730,]    4
##  [731,]    3
##  [732,]    5
##  [733,]    6
##  [734,]    5
##  [735,]    6
##  [736,]    4
##  [737,]    6
##  [738,]    6
##  [739,]    6
##  [740,]    8
##  [741,]    7
##  [742,]    6
##  [743,]    9
##  [744,]    9
##  [745,]    9
##  [746,]    8
##  [747,]    5
##  [748,]    3
##  [749,]    8
##  [750,]    7
##  [751,]    4
##  [752,]    6
##  [753,]    5
##  [754,]    7
##  [755,]    5
##  [756,]    5
##  [757,]    7
##  [758,]    5
##  [759,]    7
##  [760,]    5
##  [761,]    6
##  [762,]    9
##  [763,]    5
##  [764,]    8
##  [765,]    7
##  [766,]    7
##  [767,]    2
##  [768,]    8
##  [769,]    7
##  [770,]    4
##  [771,]    5
##  [772,]    6
##  [773,]    8
##  [774,]    7
##  [775,]    8
##  [776,]    7
##  [777,]    7
##  [778,]    5
##  [779,]    5
##  [780,]    4
##  [781,]    3
##  [782,]    7
##  [783,]    5
##  [784,]    3
##  [785,]    4
##  [786,]    3
##  [787,]    3
##  [788,]    5
##  [789,]    5
##  [790,]    8
##  [791,]    9
##  [792,]    8
##  [793,]    9
##  [794,]    8
##  [795,]    5
##  [796,]    8
##  [797,]    8
##  [798,]    6
##  [799,]    7
##  [800,]    6
##  [801,]    8
##  [802,]    9
##  [803,]    9
##  [804,]    8
##  [805,]    6
##  [806,]    8
##  [807,]    6
##  [808,]    7
##  [809,]    8
##  [810,]    8
##  [811,]    4
##  [812,]    6
##  [813,]    7
##  [814,]    9
##  [815,]    7
##  [816,]    4
##  [817,]    9
##  [818,]    6
##  [819,]    5
##  [820,]    6
##  [821,]    3
##  [822,]    5
##  [823,]    6
##  [824,]    7
##  [825,]    5
##  [826,]    6
##  [827,]    3
##  [828,]    3
##  [829,]    5
##  [830,]    5
##  [831,]    5
##  [832,]    5
##  [833,]    8
##  [834,]    4
##  [835,]    6
##  [836,]    4
##  [837,]    7
##  [838,]    1
##  [839,]    4
##  [840,]    2
##  [841,]    8
##  [842,]    8
##  [843,]    8
##  [844,]    8
##  [845,]    6
##  [846,]    3
##  [847,]    4
##  [848,]    3
##  [849,]    3
##  [850,]    6
##  [851,]    8
##  [852,]    8
##  [853,]    5
##  [854,]    7
##  [855,]    5
##  [856,]    5
##  [857,]    4
##  [858,]    8
##  [859,]    8
##  [860,]    8
##  [861,]    7
##  [862,]    7
##  [863,]    5
##  [864,]    8
##  [865,]    5
##  [866,]    7
##  [867,]    8
##  [868,]    9
##  [869,]    6
##  [870,]    8
##  [871,]    2
##  [872,]    7
##  [873,]    9
##  [874,]    4
##  [875,]    3
##  [876,]    0
##  [877,]    5
##  [878,]    4
##  [879,]    2
##  [880,]    3
##  [881,]    0
##  [882,]    9
##  [883,]    9
##  [884,]    7
##  [885,]    5
##  [886,]    5
##  [887,]    5
##  [888,]    3
##  [889,]    5
##  [890,]    6
##  [891,]    6
##  [892,]    3
##  [893,]    7
##  [894,]    6
##  [895,]    4
##  [896,]    7
##  [897,]    4
##  [898,]    8
##  [899,]    6
##  [900,]    5
##  [901,]    8
##  [902,]    6
##  [903,]    5
##  [904,]    6
##  [905,]    5
##  [906,]    7
##  [907,]    7
##  [908,]    9
##  [909,]    8
##  [910,]    8
##  [911,]    7
##  [912,]    8
##  [913,]    4
##  [914,]    9
##  [915,]    4
##  [916,]    8
##  [917,]    5
##  [918,]    5
##  [919,]    6
##  [920,]    4
##  [921,]    6
##  [922,]    5
##  [923,]    6
##  [924,]    5
##  [925,]    8
##  [926,]    6
##  [927,]    3
##  [928,]    4
##  [929,]    5
##  [930,]    5
##  [931,]    5
##  [932,]    7
##  [933,]    6
##  [934,]    6
##  [935,]    2
##  [936,]    6
##  [937,]    2
##  [938,]    7
##  [939,]    9
##  [940,]    8
##  [941,]    2
##  [942,]    4
##  [943,]    7
##  [944,]    8
##  [945,]    8
##  [946,]    7
##  [947,]    4
##  [948,]    3
##  [949,]    7
##  [950,]    4
##  [951,]    6
##  [952,]    5
##  [953,]    3
##  [954,]    5
##  [955,]    4
##  [956,]    2
##  [957,]    6
##  [958,]    6
##  [959,]    8
##  [960,]    9
##  [961,]    7
##  [962,]    7
##  [963,]    6
##  [964,]    6
##  [965,]    3
##  [966,]    3
##  [967,]    4
##  [968,]    4
##  [969,]    4
##  [970,]    5
##  [971,]    6
##  [972,]    6
##  [973,]    3
##  [974,]    6
##  [975,]    4
##  [976,]    6
##  [977,]    6
##  [978,]    9
##  [979,]    9
##  [980,]    9
##  [981,]    8
##  [982,]    8
##  [983,]    4
##  [984,]    4
##  [985,]    6
##  [986,]    9
##  [987,]    3
##  [988,]    5
##  [989,]    2
##  [990,]    5
##  [991,]    5
##  [992,]    6
##  [993,]    7
##  [994,]    7
##  [995,]    7
##  [996,]    7
##  [997,]    5
##  [998,]    1
##  [999,]    4
## [1000,]    3
## [1001,]    4
## [1002,]    5
## [1003,]    7
## [1004,]    9
## [1005,]    6
## [1006,]    9
## [1007,]    6
## [1008,]    8
## [1009,]    8
## [1010,]    5
## [1011,]    4
## [1012,]    6
## [1013,]    6
## [1014,]    7
## [1015,]    6
## [1016,]    4
## [1017,]    6
## [1018,]    4
## [1019,]    6
## [1020,]    3
## [1021,]    3
## [1022,]    7
## [1023,]    3
## [1024,]    7
## [1025,]    2
## [1026,]    7
## [1027,]    7
## [1028,]    8
## [1029,]    7
## [1030,]    7
## [1031,]    8
## [1032,]    7
## [1033,]    7
## [1034,]    4
## [1035,]    5
## [1036,]    5
## [1037,]    3
## [1038,]    6
## [1039,]    2
## [1040,]    7
## [1041,]    3
## [1042,]    7
## [1043,]    5
## [1044,]    3
## [1045,]    7
## [1046,]    6
## [1047,]    4
## [1048,]    6
## [1049,]    8
## [1050,]    9
## [1051,]    5
## [1052,]    6
## [1053,]    5
## [1054,]    8
## [1055,]    6
## [1056,]    7
## [1057,]    1
## [1058,]    5
## [1059,]    3
## [1060,]    7
## [1061,]    5
## [1062,]    8
## [1063,]    4
## [1064,]    9
## [1065,]    9
## [1066,]    4
## [1067,]    6
## [1068,]    7
## [1069,]    7
## [1070,]    8
## [1071,]    8
## [1072,]    4
## [1073,]    7
## [1074,]    5
## [1075,]    4
## [1076,]    7
## [1077,]    8
## [1078,]    3
## [1079,]    6
## [1080,]    9
## [1081,]    8
## [1082,]    6
## [1083,]    5
## [1084,]    8
## [1085,]    8
## [1086,]    8
## [1087,]    7
## [1088,]    8
## [1089,]    7
## [1090,]    5
## [1091,]    8
## [1092,]    6
## [1093,]    3
## [1094,]    6
## [1095,]    7
## [1096,]    6
## [1097,]    7
## [1098,]    8
## [1099,]    3
## [1100,]    5
## [1101,]    3
## [1102,]    8
## [1103,]    1
## [1104,]    5
## [1105,]    8
## [1106,]    7
## [1107,]    8
## [1108,]    3
## [1109,]    3
## [1110,]    7
## [1111,]    8
## [1112,]    9
## [1113,]    7
## [1114,]    7
## [1115,]    9
## [1116,]    8
## [1117,]    8
## [1118,]    6
## [1119,]    8
## [1120,]    7
## [1121,]    5
## [1122,]    7
## [1123,]    3
## [1124,]    8
## [1125,]    8
## [1126,]    5
## [1127,]    5
## [1128,]    6
## [1129,]    4
## [1130,]    3
## [1131,]    5
## [1132,]    6
## [1133,]    8
## [1134,]    5
## [1135,]    7
## [1136,]    5
## [1137,]    5
## [1138,]    6
## [1139,]    2
## [1140,]    2
## [1141,]    7
## [1142,]    4
## [1143,]    8
## [1144,]    8
## [1145,]    3
## [1146,]    9
## [1147,]    9
## [1148,]    7
## [1149,]    7
## [1150,]    6
## [1151,]    6
## [1152,]    5
## [1153,]    6
## [1154,]    9
## [1155,]    9
## [1156,]    9
## [1157,]    7
## [1158,]    2
## [1159,]    8
## [1160,]    6
## [1161,]    4
## [1162,]    7
## [1163,]    2
## [1164,]    6
## [1165,]    6
## [1166,]    3
## [1167,]    5
## [1168,]    3
## [1169,]    6
## [1170,]    8
## [1171,]    6
## [1172,]    7
## [1173,]    3
## [1174,]    2
## [1175,]    4
## [1176,]    9
## [1177,]    8
## [1178,]    9
## [1179,]    7
## [1180,]    8
## [1181,]    4
## [1182,]    6
## [1183,]    6
## [1184,]    9
## [1185,]    8
## [1186,]    8
## [1187,]    4
## [1188,]    4
## [1189,]    6
## [1190,]    6
## [1191,]    3
## [1192,]    7
## [1193,]    3
## [1194,]    5
## [1195,]    7
## [1196,]    5
## [1197,]    7
## [1198,]    8
## [1199,]    9
## [1200,]    5
## [1201,]    5
## [1202,]    9
## [1203,]    2
## [1204,]    2
## [1205,]    4
## [1206,]    3
## [1207,]    4
## [1208,]    5
## [1209,]    7
## [1210,]    8
## [1211,]    6
## [1212,]    7
## [1213,]    8
## [1214,]    8
## [1215,]    8
## [1216,]    7
## [1217,]    5
## [1218,]    7
## [1219,]    8
## [1220,]    8
## [1221,]    6
## [1222,]    8
## [1223,]    5
## [1224,]    7
## [1225,]    6
## [1226,]    6
## [1227,]    6
## [1228,]    5
## [1229,]    4
## [1230,]    6
## [1231,]    5
## [1232,]    7
## [1233,]    5
## [1234,]    1
## [1235,]    3
## [1236,]    3
## [1237,]    3
## [1238,]    2
## [1239,]    8
## [1240,]    6
## [1241,]    9
## [1242,]    9
## [1243,]    9
## [1244,]    8
## [1245,]    6
## [1246,]    9
## [1247,]    8
## [1248,]    7
## [1249,]    5
## [1250,]    5
## [1251,]    7
## [1252,]    5
## [1253,]    9
## [1254,]    7
## [1255,]    9
## [1256,]    7
## [1257,]    6
## [1258,]    5
## [1259,]    7
## [1260,]    7
## [1261,]    9
## [1262,]    6
## [1263,]    2
## [1264,]    3
## [1265,]    2
## [1266,]    3
## [1267,]    3
## [1268,]    7
## [1269,]    6
## [1270,]    6
## [1271,]    5
## [1272,]    6
## [1273,]    6
## [1274,]    4
## [1275,]    7
## [1276,]    7
## [1277,]    5
## [1278,]    4
## [1279,]    6
## [1280,]    7
## [1281,]    4
## [1282,]    8
## [1283,]    3
## [1284,]    5
## [1285,]    2
## [1286,]    4
## [1287,]    6
## [1288,]    3
## [1289,]    5
## [1290,]    8
## [1291,]    6
## [1292,]    9
## [1293,]    5
## [1294,]    7
## [1295,]    6
## [1296,]    8
## [1297,]    4
## [1298,]    8
## [1299,]    8
## [1300,]    6
## [1301,]    6
## [1302,]    5
## [1303,]    7
## [1304,]    5
## [1305,]    7
## [1306,]    9
## [1307,]    7
## [1308,]    8
## [1309,]    7
## [1310,]    8
## [1311,]    6
## [1312,]    9
## [1313,]    8
## [1314,]    5
## [1315,]    4
## [1316,]    4
## [1317,]    7
## [1318,]    4
## [1319,]    8
## [1320,]    7
## [1321,]    6
## [1322,]    7
## [1323,]    5
## [1324,]    6
## [1325,]    6
## [1326,]    7
## [1327,]    4
## [1328,]    8
## [1329,]    8
## [1330,]    6
## [1331,]    8
## [1332,]    6
## [1333,]    4
## [1334,]    5
## [1335,]    4
## [1336,]    5
## [1337,]    6
## [1338,]    3
## [1339,]    7
## [1340,]    8
## [1341,]    8
## [1342,]    9
## [1343,]    8
## [1344,]    7
## [1345,]    6
## [1346,]    6
## [1347,]    7
## [1348,]    4
## [1349,]    6
## [1350,]    4
## [1351,]    8
## [1352,]    4
## [1353,]    7
## [1354,]    6
## [1355,]    7
## [1356,]    7
## [1357,]    6
## [1358,]    3
## [1359,]    5
## [1360,]    2
## [1361,]    5
## [1362,]    3
## [1363,]    2
## [1364,]    6
## [1365,]    2
## [1366,]    6
## [1367,]    6
## [1368,]    3
## [1369,]    6
## [1370,]    5
## [1371,]    5
## [1372,]    1
## [1373,]    7
## [1374,]    4
## [1375,]    8
## [1376,]    6
## [1377,]    7
## [1378,]    2
## [1379,]    3
## [1380,]    3
## [1381,]    5
## [1382,]    6
## [1383,]    4
## [1384,]    5
## [1385,]    8
## [1386,]    7
## [1387,]    8
## [1388,]    3
## [1389,]    7
## [1390,]    6
## [1391,]    9
## [1392,]    4
## [1393,]    7
## [1394,]    4
## [1395,]    4
## [1396,]    8
## [1397,]    7
## [1398,]    6
## [1399,]    6
## [1400,]    5
## [1401,]    7
## [1402,]    3
## [1403,]    7
## [1404,]    6
## [1405,]    7
## [1406,]    5
## [1407,]    6
## [1408,]    6
## [1409,]    3
## [1410,]    9
## [1411,]    9
## [1412,]    7
## [1413,]    4
## [1414,]    6
## [1415,]    5
## [1416,]    6
## [1417,]    7
## [1418,]    8
## [1419,]    6
## [1420,]    8
## [1421,]    8
## [1422,]    8
## [1423,]    4
## [1424,]    6
## [1425,]    5
## [1426,]    6
## [1427,]    4
## [1428,]    5
## [1429,]    5
## [1430,]    6
## [1431,]    4
## [1432,]    3
## [1433,]    4
## [1434,]    4
## [1435,]    5
## [1436,]    3
## [1437,]    2
## [1438,]    7
## [1439,]    3
## [1440,]    6
## [1441,]    6
## [1442,]    4
## [1443,]    6
## [1444,]    7
## [1445,]    4
## [1446,]    6
## [1447,]    6
## [1448,]    8
## [1449,]    7
## [1450,]    2
## [1451,]    7
## [1452,]    7
## [1453,]    7
## [1454,]    3
## [1455,]    8
## [1456,]    7
## [1457,]    8
## [1458,]    5
## [1459,]    8
## [1460,]    3
## [1461,]    4
## [1462,]    6
## [1463,]    2
## [1464,]    5
## [1465,]    3
## [1466,]    4
## [1467,]    7
## [1468,]    6
## [1469,]    4
## [1470,]    6
## [1471,]    4
## [1472,]    9
## [1473,]    5
## [1474,]    7
## [1475,]    3
## [1476,]    4
## [1477,]    7
## [1478,]    5
## [1479,]    9
## [1480,]    5
## [1481,]    7
## [1482,]    8
## [1483,]    5
## [1484,]    8
## [1485,]    8
## [1486,]    9
## [1487,]    3
## [1488,]    3
## [1489,]    7
## [1490,]    6
## [1491,]    5
## [1492,]    5
## [1493,]    5
## [1494,]    6
## [1495,]    6
## [1496,]    5
## [1497,]    7
## [1498,]    6
## [1499,]    5
## [1500,]    6
## [1501,]    6
## [1502,]    7
## [1503,]    7
## [1504,]    5
## [1505,]    7
## [1506,]    8
## [1507,]    7
## [1508,]    7
## [1509,]    6
## [1510,]    4
## [1511,]    4
## [1512,]    5
## [1513,]    2
## [1514,]    6
## [1515,]    8
## [1516,]    7
## [1517,]    2
## [1518,]    5
## [1519,]    5
## [1520,]    4
## [1521,]    3
## [1522,]    1
## [1523,]    7
## [1524,]    6
## [1525,]    7
## [1526,]    6
## [1527,]    5
## [1528,]    5
## [1529,]    9
## [1530,]    8
## [1531,]    4
## [1532,]    5
## [1533,]    5
## [1534,]    6
## [1535,]    7
## [1536,]    4
## [1537,]    8
## [1538,]    5
## [1539,]    5
## [1540,]    5
## [1541,]    5
## [1542,]    4
## [1543,]    4
## [1544,]    3
## [1545,]    2
## [1546,]    5
## [1547,]    3
## [1548,]    6
## [1549,]    6
## [1550,]    7
## [1551,]    8
## [1552,]    3
## [1553,]    4
## [1554,]    5
## [1555,]    6
## [1556,]    5
## [1557,]    6
## [1558,]    5
## [1559,]    5
## [1560,]    6
## [1561,]    8
## [1562,]    5
## [1563,]    6
## [1564,]    4
## [1565,]    5
## [1566,]    6
## [1567,]    6
## [1568,]    7
## [1569,]    8
## [1570,]    5
## [1571,]    2
## [1572,]    5
## [1573,]    9
## [1574,]    9
## [1575,]    8
## [1576,]    8
## [1577,]    7
## [1578,]    5
## [1579,]    2
## [1580,]    2
## [1581,]    5
## [1582,]    4
## [1583,]    5
## [1584,]    4
## [1585,]    6
## [1586,]    7
## [1587,]    5
## [1588,]    7
## [1589,]    5
## [1590,]    6
## [1591,]    5
## [1592,]    2
## [1593,]    7
## [1594,]    5
## [1595,]    5
## [1596,]    5
## [1597,]    7
## [1598,]    7
## [1599,]    6
## [1600,]    6
## [1601,]    7
## [1602,]    6
## [1603,]    6
## [1604,]    8
## [1605,]    7
## [1606,]    7
## [1607,]    6
## [1608,]    6
## [1609,]    8
## [1610,]    6
## [1611,]    6
## [1612,]    2
## [1613,]    5
## [1614,]    5
## [1615,]    5
## [1616,]    4
## [1617,]    4
## [1618,]    7
## [1619,]    4
## [1620,]    6
## [1621,]    7
## [1622,]    7
## [1623,]    6
## [1624,]    9
## [1625,]    6
## [1626,]    2
## [1627,]    2
## [1628,]    5
## [1629,]    4
## [1630,]    8
## [1631,]    6
## [1632,]    4
## [1633,]    5
## [1634,]    8
## [1635,]    3
## [1636,]    5
## [1637,]    7
## [1638,]    3
## [1639,]    4
## [1640,]    4
## [1641,]    4
## [1642,]    7
## [1643,]    5
## [1644,]    4
## [1645,]    6
## [1646,]    7
## [1647,]    2
## [1648,]    5
## [1649,]    7
## [1650,]    8
## [1651,]    7
## [1652,]    4
## [1653,]    5
## [1654,]    8
## [1655,]    8
## [1656,]    5
## [1657,]    6
## [1658,]    6
## [1659,]    6
## [1660,]    4
## [1661,]    5
## [1662,]    1
## [1663,]    4
## [1664,]    3
## [1665,]    2
## [1666,]    6
## [1667,]    3
## [1668,]    2
## [1669,]    8
## [1670,]    8
## [1671,]    4
## [1672,]    8
## [1673,]    5
## [1674,]    6
## [1675,]    5
## [1676,]    5
## [1677,]    7
## [1678,]    6
## [1679,]    2
## [1680,]    1
## [1681,]    6
## [1682,]    5
## [1683,]    6
## [1684,]    6
## [1685,]    7
## [1686,]    7
## [1687,]    4
## [1688,]    9
## [1689,]    4
## [1690,]    5
## [1691,]    5
## [1692,]    5
## [1693,]    2
## [1694,]    5
## [1695,]    8
## [1696,]    6
## [1697,]    8
## [1698,]    4
## [1699,]    8
## [1700,]    8
## [1701,]    5
## [1702,]    5
## [1703,]    9
## [1704,]    7
## [1705,]    7
## [1706,]    2
## [1707,]    8
## [1708,]    5
## [1709,]    7
## [1710,]    6
## [1711,]    8
## [1712,]    5
## [1713,]    5
## [1714,]    4
## [1715,]    8
## [1716,]    5
## [1717,]    3
## [1718,]    4
## [1719,]    6
## [1720,]    8
## [1721,]    7
## [1722,]    8
## [1723,]    7
## [1724,]    5
## [1725,]    5
## [1726,]    5
## [1727,]    7
## [1728,]    8
## [1729,]    7
## [1730,]    5
## [1731,]    5
## [1732,]    1
## [1733,]    5
## [1734,]    5
## [1735,]    5
## [1736,]    8
## [1737,]    6
## [1738,]    5
## [1739,]    5
## [1740,]    4
## [1741,]    9
## [1742,]    6
## [1743,]    7
## [1744,]    5
## [1745,]    7
## [1746,]    7
## [1747,]    8
## [1748,]    7
## [1749,]    9
## [1750,]    8
## [1751,]    5
## [1752,]    5
## [1753,]    7
## [1754,]    5
## [1755,]    5
## [1756,]    7
## [1757,]    9
## [1758,]    2
## [1759,]    2
## [1760,]    3
## [1761,]    2
## [1762,]    7
## [1763,]    6
## [1764,]    6
## [1765,]    7
## [1766,]    6
## [1767,]    7
## [1768,]    2
## [1769,]    4
## [1770,]    6
## [1771,]    7
## [1772,]    4
## [1773,]    4
## [1774,]    2
## [1775,]    5
## [1776,]    6
## [1777,]    6
## [1778,]    7
## [1779,]    7
## [1780,]    5
## [1781,]    8
## [1782,]    8
## [1783,]    8
## [1784,]    9
## [1785,]    9
## [1786,]    7
## [1787,]    6
## [1788,]    5
## [1789,]    5
## [1790,]    6
## [1791,]    5
## [1792,]    6
## [1793,]    5
## [1794,]    5
## [1795,]    8
## [1796,]    6
## [1797,]    9
## [1798,]    6
## [1799,]    4
## [1800,]    5
## [1801,]    3
## [1802,]    6
## [1803,]    8
## [1804,]    7
## [1805,]    8
## [1806,]    7
## [1807,]    5
## [1808,]    9
## [1809,]    6
## [1810,]    7
## [1811,]    8
## [1812,]    7
## [1813,]    3
## [1814,]    5
## [1815,]    6
## [1816,]    5
## [1817,]    2
## [1818,]    8
## [1819,]    2
## [1820,]    3
## [1821,]    6
## [1822,]    3
## [1823,]    2
## [1824,]    7
## [1825,]    4
## [1826,]    8
## [1827,]    6
## [1828,]    6
## [1829,]    7
## [1830,]    7
## [1831,]    8
## [1832,]    7
## [1833,]    7
## [1834,]    2
## [1835,]    8
## [1836,]    3
## [1837,]    2
## [1838,]    5
## [1839,]    8
## [1840,]    5
## [1841,]    2
## [1842,]    5
## [1843,]    7
## [1844,]    2
## [1845,]    7
## [1846,]    4
## [1847,]    8
## [1848,]    7
## [1849,]    8
## [1850,]    5
## [1851,]    7
## [1852,]    7
## [1853,]    6
## [1854,]    2
## [1855,]    7
## [1856,]    7
## [1857,]    4
## [1858,]    3
## [1859,]    6
## [1860,]    7
## [1861,]    7
## [1862,]    5
## [1863,]    2
## [1864,]    4
## [1865,]    7
## [1866,]    4
## [1867,]    4
## [1868,]    8
## [1869,]    8
## [1870,]    8
## [1871,]    7
## [1872,]    7
## [1873,]    5
## [1874,]    9
## [1875,]    8
## [1876,]    7
## [1877,]    8
## [1878,]    3
## [1879,]    5
## [1880,]    4
## [1881,]    8
## [1882,]    9
## [1883,]    9
## [1884,]    7
## [1885,]    3
## [1886,]    6
## [1887,]    4
## [1888,]    7
## [1889,]    8
## [1890,]    5
## [1891,]    5
## [1892,]    7
## [1893,]    9
## [1894,]    5
## [1895,]    7
## [1896,]    6
## [1897,]    2
## [1898,]    7
## [1899,]    6
## [1900,]    6
## [1901,]    5
## [1902,]    5
## [1903,]    7
## [1904,]    8
## [1905,]    5
## [1906,]    7
## [1907,]    4
## [1908,]    6
## [1909,]    7
## [1910,]    9
## [1911,]    8
## [1912,]    7
## [1913,]    8
## [1914,]    5
## [1915,]    6
## [1916,]    7
## [1917,]    7
## [1918,]    5
## [1919,]    5
## [1920,]    8
## [1921,]    4
## [1922,]    6
## [1923,]    5
## [1924,]    4
## [1925,]    9
## [1926,]    7
## [1927,]    4
## [1928,]    4
## [1929,]    9
## [1930,]    5
## [1931,]    5
## [1932,]    6
## [1933,]    8
## [1934,]    5
## [1935,]    6
## [1936,]    7
## [1937,]    2
## [1938,]    2
## [1939,]    3
## [1940,]    4
## [1941,]    3
## [1942,]    3
## [1943,]    4
## [1944,]    3
## [1945,]    8
## [1946,]    5
## [1947,]    4
## [1948,]    3
## [1949,]    5
## [1950,]    5
## [1951,]    8
## [1952,]    8
## [1953,]    5
## [1954,]    4
## [1955,]    3
## [1956,]    7
## [1957,]    8
## [1958,]    7
## [1959,]    5
## [1960,]    9
## [1961,]    8
## [1962,]    1
## [1963,]    8
## [1964,]    4
## [1965,]    7
## [1966,]    6
## [1967,]    6
## [1968,]    6
## [1969,]    6
## [1970,]    7
## [1971,]    5
## [1972,]    6
## [1973,]    3
## [1974,]    5
## [1975,]    7
## [1976,]    7
## [1977,]    5
## [1978,]    5
## [1979,]    5
## [1980,]    5
## [1981,]    6
## [1982,]    3
## [1983,]    5
## [1984,]    5
## [1985,]    7
## [1986,]    8
## [1987,]    7
## [1988,]    8
## [1989,]    7
## [1990,]    9
## [1991,]    5
## [1992,]    5
## [1993,]    6
## [1994,]    5
## [1995,]    1
## [1996,]    3
## [1997,]    2
## [1998,]    2
## [1999,]    4
## [2000,]    9
## [2001,]    8
## [2002,]    4
## [2003,]    6
## [2004,]    6
## [2005,]    5
## [2006,]    6
## [2007,]    8
## [2008,]    9
## [2009,]    7
## [2010,]    6
## [2011,]    3
## [2012,]    6
## [2013,]    7
## [2014,]    8
## [2015,]    6
## [2016,]    9
## [2017,]    6
## [2018,]    5
## [2019,]    6
## [2020,]    5
## [2021,]    8
## [2022,]    2
## [2023,]    8
## [2024,]    9
## [2025,]    5
## [2026,]    7
## [2027,]    6
## [2028,]    2
## [2029,]    7
## [2030,]    6
## [2031,]    5
## [2032,]    9
## [2033,]    6
## [2034,]    5
## [2035,]    5
## [2036,]    4
## [2037,]    5
## [2038,]    5
## [2039,]    6
## [2040,]    5
## [2041,]    5
## [2042,]    3
## [2043,]    6
## [2044,]    7
## [2045,]    4
## [2046,]    7
## [2047,]    7
## [2048,]    5
## [2049,]    2
## [2050,]    2
## [2051,]    5
## [2052,]    3
## [2053,]    5
## [2054,]    2
## [2055,]    3
## [2056,]    8
## [2057,]    6
## [2058,]    5
## [2059,]    5
## [2060,]    4
## [2061,]    8
## [2062,]    3
## [2063,]    2
## [2064,]    8
## [2065,]    8
## [2066,]    8
## [2067,]    6
## [2068,]    6
## [2069,]    5
## [2070,]    9
## [2071,]    8
## [2072,]    8
## [2073,]    5
## [2074,]    4
## [2075,]    4
## [2076,]    2
## [2077,]    4
## [2078,]    3
## [2079,]    3
## [2080,]    4
## [2081,]    5
## [2082,]    3
## [2083,]    5
## [2084,]    6
## [2085,]    7
## [2086,]    5
## [2087,]    6
## [2088,]    6
## [2089,]    7
## [2090,]    5
## [2091,]    6
## [2092,]    8
## [2093,]    6
## [2094,]    4
## [2095,]    4
## [2096,]    5
## [2097,]    5
## [2098,]    9
## [2099,]    9
## [2100,]    7
## [2101,]    4
## [2102,]    7
## [2103,]    5
## [2104,]    3
## [2105,]    5
## [2106,]    7
## [2107,]    7
## [2108,]    8
## [2109,]    8
## [2110,]    7
## [2111,]    7
## [2112,]    7
## [2113,]    7
## [2114,]    8
## [2115,]    7
## [2116,]    6
## [2117,]    8
## [2118,]    8
## [2119,]    6
## [2120,]    4
## [2121,]    5
## [2122,]    4
## [2123,]    2
## [2124,]    1
## [2125,]    5
## [2126,]    4
## [2127,]    6
## [2128,]    7
## [2129,]    6
## [2130,]    3
## [2131,]    4
## [2132,]    3
## [2133,]    5
## [2134,]    6
## [2135,]    5
## [2136,]    5
## [2137,]    7
## [2138,]    7
## [2139,]    7
## [2140,]    6
## [2141,]    7
## [2142,]    6
## [2143,]    7
## [2144,]    4
## [2145,]    7
## [2146,]    7
## [2147,]    3
## [2148,]    5
## [2149,]    6
## [2150,]    5
## [2151,]    8
## [2152,]    6
## [2153,]    4
## [2154,]    7
## [2155,]    5
## [2156,]    8
## [2157,]    6
## [2158,]    5
## [2159,]    6
## [2160,]    8
## [2161,]    4
## [2162,]    6
## [2163,]    6
## [2164,]    9
## [2165,]    5
## [2166,]    4
## [2167,]    5
## [2168,]    4
## [2169,]    7
## [2170,]    7
## [2171,]    6
## [2172,]    2
## [2173,]    8
## [2174,]    6
## [2175,]    7
## [2176,]    4
## [2177,]    5
## [2178,]    6
## [2179,]    7
## [2180,]    4
## [2181,]    7
## [2182,]    5
## [2183,]    3
## [2184,]    5
## [2185,]    3
## [2186,]    3
## [2187,]    5
## [2188,]    8
## [2189,]    9
## [2190,]    9
## [2191,]    7
## [2192,]    6
## [2193,]    7
## [2194,]    5
## [2195,]    8
## [2196,]    6
## [2197,]    3
## [2198,]    9
## [2199,]    9
## [2200,]    8
## [2201,]    7
## [2202,]    5
## [2203,]    6
## [2204,]    6
## [2205,]    4
## [2206,]    5
## [2207,]    5
## [2208,]    6
## [2209,]    5
## [2210,]    6
## [2211,]    4
## [2212,]    8
## [2213,]    4
## [2214,]    3
## [2215,]    4
## [2216,]    2
## [2217,]    5
## [2218,]    7
## [2219,]    4
## [2220,]    6
## [2221,]    6
## [2222,]    4
## [2223,]    5
## [2224,]    3
## [2225,]    4
## [2226,]    5
## [2227,]    5
## [2228,]    6
## [2229,]    4
## [2230,]    4
## [2231,]    5
## [2232,]    7
## [2233,]    7
## [2234,]    6
## [2235,]    5
## [2236,]    5
## [2237,]    4
## [2238,]    6
## [2239,]    5
## [2240,]    8
## [2241,]    9
## [2242,]    7
## [2243,]    6
## [2244,]    7
## [2245,]    4
## [2246,]    6
## [2247,]    6
## [2248,]    8
## [2249,]    5
## [2250,]    9
## [2251,]    7
## [2252,]    5
## [2253,]    8
## [2254,]    4
## [2255,]    4
## [2256,]    3
## [2257,]    5
## [2258,]    8
## [2259,]    5
## [2260,]    8
## [2261,]    5
## [2262,]    6
## [2263,]    7
## [2264,]    6
## [2265,]    7
## [2266,]    8
## [2267,]    5
## [2268,]    4
## [2269,]    5
## [2270,]    4
## [2271,]    8
## [2272,]    6
## [2273,]    7
## [2274,]    4
## [2275,]    7
## [2276,]    9
## [2277,]    6
## [2278,]    6
## [2279,]    4
## [2280,]    4
## [2281,]    5
## [2282,]    7
## [2283,]    9
## [2284,]    4
## [2285,]    5
## [2286,]    6
## [2287,]    9
## [2288,]    7
## [2289,]    8
## [2290,]    5
## [2291,]    7
## [2292,]    8
## [2293,]    6
## [2294,]    8
## [2295,]    8
## [2296,]    7
## [2297,]    5
## [2298,]    7
## [2299,]    6
## [2300,]    8
## [2301,]    5
## [2302,]    6
## [2303,]    7
## [2304,]    5
## [2305,]    5
## [2306,]    5
## [2307,]    4
## [2308,]    3
## [2309,]    6
## [2310,]    6
## [2311,]    5
## [2312,]    7
## [2313,]    5
## [2314,]    5
## [2315,]    2
## [2316,]    4
## [2317,]    3
## [2318,]    5
## [2319,]    4
## [2320,]    6
## [2321,]    5
## [2322,]    5
## [2323,]    5
## [2324,]    6
## [2325,]    7
## [2326,]    7
## [2327,]    6
## [2328,]    6
## [2329,]    9
## [2330,]    5
## [2331,]    7
## [2332,]    4
## [2333,]    6
## [2334,]    6
## [2335,]    8
## [2336,]    6
## [2337,]    6
## [2338,]    6
## [2339,]    6
## [2340,]    4
## [2341,]    8
## [2342,]    9
## [2343,]    9
## [2344,]    9
## [2345,]    6
## [2346,]    3
## [2347,]    8
## [2348,]    8
## [2349,]    5
## [2350,]    6
## [2351,]    4
## [2352,]    8
## [2353,]    5
## [2354,]    4
## [2355,]    4
## [2356,]    7
## [2357,]    7
## [2358,]    5
## [2359,]    5
## [2360,]    3
## [2361,]    6
## [2362,]    8
## [2363,]    8
## [2364,]    7
## [2365,]    4
## [2366,]    1
## [2367,]    3
## [2368,]    5
## [2369,]    5
## [2370,]    6
## [2371,]    8
## [2372,]    9
## [2373,]    8
## [2374,]    1
## [2375,]    6
## [2376,]    6
## [2377,]    7
## [2378,]    6
## [2379,]    5
## [2380,]    4
## [2381,]    6
## [2382,]    8
## [2383,]    9
## [2384,]    8
## [2385,]    5
## [2386,]    6
## [2387,]    3
## [2388,]    6
## [2389,]    1
## [2390,]    4
## [2391,]    8
## [2392,]    7
## [2393,]    4
## [2394,]    3
## [2395,]    9
## [2396,]    7
## [2397,]    5
## [2398,]    7
## [2399,]    4
## [2400,]    5
## [2401,]    8
## [2402,]    8
## [2403,]    8
## [2404,]    7
## [2405,]    8
## [2406,]    5
## [2407,]    8
## [2408,]    6
## [2409,]    5
## [2410,]    3
## [2411,]    5
## [2412,]    7
## [2413,]    7
## [2414,]    3
## [2415,]    6
## [2416,]    6
## [2417,]    7
## [2418,]    6
## [2419,]    8
## [2420,]    7
## [2421,]    4
## [2422,]    4
## [2423,]    6
## [2424,]    4
## [2425,]    6
## [2426,]    3
## [2427,]    5
## [2428,]    6
## [2429,]    7
## [2430,]    5
## [2431,]    7
## [2432,]    3
## [2433,]    7
## [2434,]    4
## [2435,]    6
## [2436,]    8
## [2437,]    5
## [2438,]    5
## [2439,]    5
## [2440,]    5
## [2441,]    3
## [2442,]    4
## [2443,]    4
## [2444,]    3
## [2445,]    4
## [2446,]    3
## [2447,]    4
## [2448,]    6
## [2449,]    5
## [2450,]    4
## [2451,]    4
## [2452,]    7
## [2453,]    3
## [2454,]    4
## [2455,]    9
## [2456,]    6
## [2457,]    9
## [2458,]    4
## [2459,]    7
## [2460,]    5
## [2461,]    5
## [2462,]    7
## [2463,]    7
## [2464,]    3
## [2465,]    3
## [2466,]    4
## [2467,]    4
## [2468,]    8
## [2469,]    8
## [2470,]    3
## [2471,]    6
## [2472,]    5
## [2473,]    6
## [2474,]    5
## [2475,]    7
## [2476,]    8
## [2477,]    8
## [2478,]    7
## [2479,]    7
## [2480,]    5
## [2481,]    6
## [2482,]    8
## [2483,]    7
## [2484,]    6
## [2485,]    4
## [2486,]    5
## [2487,]    5
## [2488,]    5
## [2489,]    6
## [2490,]    7
## [2491,]    3
## [2492,]    7
## [2493,]    4
## [2494,]    6
## [2495,]    7
## [2496,]    7
## [2497,]    3
## [2498,]    6
## [2499,]    7
## [2500,]    6
## [2501,]    7
## [2502,]    6
## [2503,]    6
## [2504,]    7
## [2505,]    8
## [2506,]    6
## [2507,]    7
## [2508,]    9
## [2509,]    5
## [2510,]    7
## [2511,]    8
## [2512,]    3
## [2513,]    6
## [2514,]    9
## [2515,]    8
## [2516,]    7
## [2517,]    6
## [2518,]    8
## [2519,]    8
## [2520,]    7
## [2521,]    6
## [2522,]    6
## [2523,]    6
## [2524,]    8
## [2525,]    6
## [2526,]    7
## [2527,]    4
## [2528,]    7
## [2529,]    6
## [2530,]    8
## [2531,]    7
## [2532,]    7
## [2533,]    9
## [2534,]    7
## [2535,]    3
## [2536,]    6
## [2537,]    6
## [2538,]    4
## [2539,]    8
## [2540,]    8
## [2541,]    6
## [2542,]    5
## [2543,]    2
## [2544,]    6
## [2545,]    3
## [2546,]    3
## [2547,]    5
## [2548,]    8
## [2549,]    5
## [2550,]    4
## [2551,]    5
## [2552,]    6
## [2553,]    6
## [2554,]    5
## [2555,]    4
## [2556,]    8
## [2557,]    6
## [2558,]    8
## [2559,]    8
## [2560,]    3
## [2561,]    8
## [2562,]    8
## [2563,]    4
## [2564,]    8
## [2565,]    7
## [2566,]    7
## [2567,]    8
## [2568,]    8
## [2569,]    6
## [2570,]    6
## [2571,]    8
## [2572,]    5
## [2573,]    7
## [2574,]    6
## [2575,]    3
## [2576,]    7
## [2577,]    5
## [2578,]    7
## [2579,]    9
## [2580,]    3
## [2581,]    6
## [2582,]    4
## [2583,]    9
## [2584,]    4
## [2585,]    6
## [2586,]    7
## [2587,]    8
## [2588,]    7
## [2589,]    8
## [2590,]    6
## [2591,]    8
## [2592,]    4
## [2593,]    6
## [2594,]    3
## [2595,]    5
## [2596,]    6
## [2597,]    5
## [2598,]    6
## [2599,]    5
## [2600,]    2
## [2601,]    4
## [2602,]    5
## [2603,]    4
## [2604,]    6
## [2605,]    5
## [2606,]    5
## [2607,]    5
## [2608,]    5
## [2609,]    9
## [2610,]    4
## [2611,]    7
## [2612,]    5
## [2613,]    1
## [2614,]    7
## [2615,]    6
## [2616,]    5
## [2617,]    3
## [2618,]    5
## [2619,]    6
## [2620,]    7
## [2621,]    3
## [2622,]    5
## [2623,]    5
## [2624,]    4
## [2625,]    5
## [2626,]    6
## [2627,]    5
## [2628,]    6
## [2629,]    8
## [2630,]    8
## [2631,]    5
## [2632,]    4
## [2633,]    4
## [2634,]    4
## [2635,]    4
## [2636,]    6
## [2637,]    7
## [2638,]    4
## [2639,]    7
## [2640,]    6
## [2641,]    7
## [2642,]    8
## [2643,]    8
## [2644,]    8
## [2645,]    7
## [2646,]    7
## [2647,]    4
## [2648,]    4
## [2649,]    9
## [2650,]    7
## [2651,]    3
## [2652,]    2
## [2653,]    2
## [2654,]    3
## [2655,]    4
## [2656,]    6
## [2657,]    5
## [2658,]    2
## [2659,]    2
## [2660,]    7
## [2661,]    5
## [2662,]    7
## [2663,]    5
## [2664,]    6
## [2665,]    6
## [2666,]    6
## [2667,]    8
## [2668,]    7
## [2669,]    8
## [2670,]    8
## [2671,]    4
## [2672,]    8
## [2673,]    4
## [2674,]    5
## [2675,]    7
## [2676,]    7
## [2677,]    6
## [2678,]    5
## [2679,]    7
## [2680,]    6
## [2681,]    6
## [2682,]    8
## [2683,]    4
## [2684,]    7
## [2685,]    7
## [2686,]    6
## [2687,]    6
## [2688,]    5
## [2689,]    9
## [2690,]    7
## [2691,]    5
## [2692,]    7
## [2693,]    7
## [2694,]    8
## [2695,]    5
## [2696,]    6
## [2697,]    3
## [2698,]    1
## [2699,]    5
## [2700,]    4
## [2701,]    5
## [2702,]    2
## [2703,]    4
## [2704,]    5
## [2705,]    2
## [2706,]    1
## [2707,]    3
## [2708,]    6
## [2709,]    4
## [2710,]    8
## [2711,]    5
## [2712,]    4
## [2713,]    5
## [2714,]    5
## [2715,]    8
## [2716,]    6
## [2717,]    6
## [2718,]    5
## [2719,]    4
## [2720,]    5
## [2721,]    7
## [2722,]    3
## [2723,]    3
## [2724,]    4
## [2725,]    8
## [2726,]    8
## [2727,]    6
## [2728,]    7
## [2729,]    5
## [2730,]    5
## [2731,]    6
## [2732,]    6
## [2733,]    7
## [2734,]    7
## [2735,]    5
## [2736,]    7
## [2737,]    5
## [2738,]    5
## [2739,]    7
## [2740,]    6
## [2741,]    4
## [2742,]    4
## [2743,]    3
## [2744,]    4
## [2745,]    2
## [2746,]    8
## [2747,]    6
## [2748,]    6
## [2749,]    8
## [2750,]    6
## [2751,]    8
## [2752,]    5
## [2753,]    3
## [2754,]    5
## [2755,]    5
## [2756,]    7
## [2757,]    9
## [2758,]    5
## [2759,]    7
## [2760,]    5
## [2761,]    8
## [2762,]    9
## [2763,]    8
## [2764,]    7
## [2765,]    8
## [2766,]    6
## [2767,]    5
## [2768,]    6
## [2769,]    4
## [2770,]    4
## [2771,]    6
## [2772,]    5
## [2773,]    7
## [2774,]    5
## [2775,]    3
## [2776,]    9
## [2777,]    8
## [2778,]    5
## [2779,]    8
## [2780,]    6
## [2781,]    8
## [2782,]    5
## [2783,]    6
## [2784,]    6
## [2785,]    3
## [2786,]    6
## [2787,]    5
## [2788,]    6
## [2789,]    4
## [2790,]    8
## [2791,]    4
## [2792,]    3
## [2793,]    5
## [2794,]    1
## [2795,]    4
## [2796,]    3
## [2797,]    7
## [2798,]    6
## [2799,]    6
## [2800,]    7
## [2801,]    5
## [2802,]    6
## [2803,]    5
## [2804,]    3
## [2805,]    5
## [2806,]    7
## [2807,]    8
## [2808,]    7
## [2809,]    6
## [2810,]    5
## [2811,]    2
## [2812,]    5
## [2813,]    6
## [2814,]    6
## [2815,]    9
## [2816,]    7
## [2817,]    8
## [2818,]    7
## [2819,]    5
## [2820,]    7
## [2821,]    8
## [2822,]    8
## [2823,]    7
## [2824,]    5
## [2825,]    6
## [2826,]    6
## [2827,]    1
## [2828,]    4
## [2829,]    3
## [2830,]    4
## [2831,]    7
## [2832,]    7
## [2833,]    7
## [2834,]    6
## [2835,]    7
## [2836,]    8
## [2837,]    7
## [2838,]    6
## [2839,]    6
## [2840,]    5
## [2841,]    5
## [2842,]    5
## [2843,]    6
## [2844,]    6
## [2845,]    6
## [2846,]    5
## [2847,]    4
## [2848,]    4
## [2849,]    5
## [2850,]    6
## [2851,]    5
## [2852,]    6
## [2853,]    5
## [2854,]    7
## [2855,]    4
## [2856,]    5
## [2857,]    4
## [2858,]    2
## [2859,]    3
## [2860,]    6
## [2861,]    5
## [2862,]    3
## [2863,]    5
## [2864,]    3
## [2865,]    6
## [2866,]    3
## [2867,]    4
## [2868,]    8
## [2869,]    7
## [2870,]    3
## [2871,]    5
## [2872,]    5
## [2873,]    7
## [2874,]    6
## [2875,]    2
## [2876,]    2
## [2877,]    3
## [2878,]    3
## [2879,]    8
## [2880,]    8
## [2881,]    6
## [2882,]    4
## [2883,]    7
## [2884,]    8
## [2885,]    7
## [2886,]    7
## [2887,]    6
## [2888,]    4
## [2889,]    6
## [2890,]    6
## [2891,]    8
## [2892,]    3
## [2893,]    3
## [2894,]    8
## [2895,]    6
## [2896,]    6
## [2897,]    5
## [2898,]    3
## [2899,]    5
## [2900,]    4
## [2901,]    8
## [2902,]    7
## [2903,]    5
## [2904,]    6
## [2905,]    9
## [2906,]    7
## [2907,]    7
## [2908,]    5
## [2909,]    5
## [2910,]    3
## [2911,]    3
## [2912,]    0
## [2913,]    5
## [2914,]    1
## [2915,]    3
## [2916,]    8
## [2917,]    8
## [2918,]    5
## [2919,]    6
## [2920,]    5
## [2921,]    7
## [2922,]    7
## [2923,]    7
## [2924,]    8
## [2925,]    6
## [2926,]    5
## [2927,]    3
## [2928,]    5
## [2929,]    8
## [2930,]    4
## [2931,]    6
## [2932,]    7
## [2933,]    4
## [2934,]    5
## [2935,]    6
## [2936,]    3
## [2937,]    5
## [2938,]    3
## [2939,]    6
## [2940,]    7
## [2941,]    8
## [2942,]    6
## [2943,]    4
## [2944,]    5
## [2945,]    6
## [2946,]    5
## [2947,]    3
## [2948,]    8
## [2949,]    6
## [2950,]    6
## [2951,]    6
## [2952,]    5
## [2953,]    7
## [2954,]    7
## [2955,]    8
## [2956,]    6
## [2957,]    7
## [2958,]    3
## [2959,]    5
## [2960,]    3
## [2961,]    6
## [2962,]    5
## [2963,]    3
## [2964,]    7
## [2965,]    8
## [2966,]    7
## [2967,]    9
## [2968,]    9
## [2969,]    8
## [2970,]    8
## [2971,]    5
## [2972,]    4
## [2973,]    2
## [2974,]    7
## [2975,]    7
## [2976,]    8
## [2977,]    4
## [2978,]    7
## [2979,]    7
## [2980,]    1
## [2981,]    7
## [2982,]    2
## [2983,]    4
## [2984,]    2
## [2985,]    8
## [2986,]    5
## [2987,]    6
## [2988,]    9
## [2989,]    4
## [2990,]    7
## [2991,]    5
## [2992,]    5
## [2993,]    2
## [2994,]    6
## [2995,]    4
## [2996,]    1
## [2997,]    3
## [2998,]    6
## [2999,]    6
## [3000,]    4
## [3001,]    4
## [3002,]    5
## [3003,]    8
## [3004,]    7
## [3005,]    8
## [3006,]    3
## [3007,]    4
## [3008,]    4
## [3009,]    5
## [3010,]    6
## [3011,]    3
## [3012,]    3
## [3013,]    9
## [3014,]    6
## [3015,]    7
## [3016,]    7
## [3017,]    5
## [3018,]    7
## [3019,]    5
## [3020,]    6
## [3021,]    3
## [3022,]    7
## [3023,]    5
## [3024,]    3
## [3025,]    4
## [3026,]    4
## [3027,]    3
## [3028,]    7
## [3029,]    7
## [3030,]    7
## [3031,]    8
## [3032,]    7
## [3033,]    8
## [3034,]    5
## [3035,]    2
## [3036,]    3
## [3037,]    6
## [3038,]    8
## [3039,]    6
## [3040,]    4
## [3041,]    7
## [3042,]    3
## [3043,]    5
## [3044,]    6
## [3045,]    5
## [3046,]    7
## [3047,]    6
## [3048,]    8
## [3049,]    6
## [3050,]    7
## [3051,]    6
## [3052,]    5
## [3053,]    5
## [3054,]    4
## [3055,]    4
## [3056,]    7
## [3057,]    6
## [3058,]    7
## [3059,]    1
## [3060,]    1
## [3061,]    6
## [3062,]    6
## [3063,]    6
## [3064,]    3
## [3065,]    6
## [3066,]    6
## [3067,]    7
## [3068,]    7
## [3069,]    3
## [3070,]    6
## [3071,]    7
## [3072,]    8
## [3073,]    7
## [3074,]    7
## [3075,]    9
## [3076,]    7
## [3077,]    6
## [3078,]    6
## [3079,]    9
## [3080,]    7
## [3081,]    7
## [3082,]    4
## [3083,]    4
## [3084,]    5
## [3085,]    2
## [3086,]    7
## [3087,]    4
## [3088,]    6
## [3089,]    8
## [3090,]    9
## [3091,]    6
## [3092,]    5
## [3093,]    7
## [3094,]    4
## [3095,]    5
## [3096,]    6
## [3097,]    4
## [3098,]    8
## [3099,]    7
## [3100,]    3
## [3101,]    3
## [3102,]    1
## [3103,]    9
## [3104,]    7
## [3105,]    7
## [3106,]    4
## [3107,]    7
## [3108,]    4
## [3109,]    7
## [3110,]    5
## [3111,]    5
## [3112,]    7
## [3113,]    6
## [3114,]    9
## [3115,]    5
## [3116,]    5
## [3117,]    5
## [3118,]    6
## [3119,]    7
## [3120,]    8
## [3121,]    4
## [3122,]    9
## [3123,]    6
## [3124,]    3
## [3125,]    2
## [3126,]    6
## [3127,]    6
## [3128,]    7
## [3129,]    4
## [3130,]    4
## [3131,]    5
## [3132,]    3
## [3133,]    6
## [3134,]    8
## [3135,]    4
## [3136,]    7
## [3137,]    5
## [3138,]    5
## [3139,]    5
## [3140,]    6
## [3141,]    5
## [3142,]    7
## [3143,]    4
## [3144,]    1
## [3145,]    4
## [3146,]    4
## [3147,]    5
## [3148,]    5
## [3149,]    8
## [3150,]    8
## [3151,]    5
## [3152,]    5
## [3153,]    5
## [3154,]    3
## [3155,]    6
## [3156,]    5
## [3157,]    4
## [3158,]    3
## [3159,]    4
## [3160,]    3
## [3161,]    8
## [3162,]    4
## [3163,]    9
## [3164,]    3
## [3165,]    6
## [3166,]    9
## [3167,]    6
## [3168,]    5
## [3169,]    6
## [3170,]    8
## [3171,]    7
## [3172,]    5
## [3173,]    8
## [3174,]    8
## [3175,]    4
## [3176,]    7
## [3177,]    6
## [3178,]    2
## [3179,]    2
## [3180,]    6
## [3181,]    6
## [3182,]    8
## [3183,]    7
## [3184,]    6
## [3185,]    7
## [3186,]    5
## [3187,]    8
## [3188,]    7
## [3189,]    7
## [3190,]    7
## [3191,]    4
## [3192,]    4
## [3193,]    5
## [3194,]    3
## [3195,]    4
## [3196,]    7
## [3197,]    6
## [3198,]    5
## [3199,]    6
## [3200,]    8
## [3201,]    4
## [3202,]    7
## [3203,]    9
## [3204,]    7
## [3205,]    6
## [3206,]    7
## [3207,]    5
## [3208,]    7
## [3209,]    5
## [3210,]    2
## [3211,]    6
## [3212,]    3
## [3213,]    6
## [3214,]    5
## [3215,]    7
## [3216,]    6
## [3217,]    6
## [3218,]    5
## [3219,]    7
## [3220,]    4
## [3221,]    3
## [3222,]    6
## [3223,]    5
## [3224,]    2
## [3225,]    5
## [3226,]    7
## [3227,]    6
## [3228,]    2
## [3229,]    7
## [3230,]    4
## [3231,]    3
## [3232,]    2
## [3233,]    3
## [3234,]    8
## [3235,]    6
## [3236,]    9
## [3237,]    3
## [3238,]    8
## [3239,]    6
## [3240,]    4
## [3241,]    4
## [3242,]    6
## [3243,]    6
## [3244,]    7
## [3245,]    8
## [3246,]    2
## [3247,]    9
## [3248,]    7
## [3249,]    6
## [3250,]    6
## [3251,]    4
## [3252,]    8
## [3253,]    8
## [3254,]    7
## [3255,]    4
## [3256,]    7
## [3257,]    6
## [3258,]    7
## [3259,]    5
## [3260,]    8
## [3261,]    6
## [3262,]    6
## [3263,]    6
## [3264,]    9
## [3265,]    6
## [3266,]    6
## [3267,]    7
## [3268,]    7
## [3269,]    5
## [3270,]    4
## [3271,]    4
## [3272,]    5
## [3273,]    5
## [3274,]    6
## [3275,]    9
## [3276,]    8
## [3277,]    7
## [3278,]    3
## [3279,]    1
## [3280,]    9
## [3281,]    8
## [3282,]    8
## [3283,]    8
## [3284,]    7
## [3285,]    7
## [3286,]    6
## [3287,]    2
## [3288,]    5
## [3289,]    5
## [3290,]    5
## [3291,]    7
## [3292,]    7
## [3293,]    8
## [3294,]    4
## [3295,]    5
## [3296,]    5
## [3297,]    5
## [3298,]    6
## [3299,]    6
## [3300,]    6
## [3301,]    7
## [3302,]    8
## [3303,]    6
## [3304,]    8
## [3305,]    5
## [3306,]    4
## [3307,]    4
## [3308,]    7
## [3309,]    7
## [3310,]    4
## [3311,]    4
## [3312,]    5
## [3313,]    4
## [3314,]    4
## [3315,]    4
## [3316,]    4
## [3317,]    7
## [3318,]    6
## [3319,]    5
## [3320,]    7
## [3321,]    6
## [3322,]    3
## [3323,]    8
## [3324,]    4
## [3325,]    7
## [3326,]    5
## [3327,]    5
## [3328,]    5
## [3329,]    5
## [3330,]    5
## [3331,]    5
## [3332,]    7
## [3333,]    6
## [3334,]    4
## [3335,]    8
## [3336,]    7
## [3337,]    7
## [3338,]    6
## [3339,]    7
## [3340,]    6
## [3341,]    8
## [3342,]    5
## [3343,]    6
## [3344,]    7
## [3345,]    5
## [3346,]    4
## [3347,]    1
## [3348,]    8
## [3349,]    8
## [3350,]    5
## [3351,]    4
## [3352,]    5
## [3353,]    5
## [3354,]    3
## [3355,]    1
## [3356,]    6
## [3357,]    5
## [3358,]    6
## [3359,]    6
## [3360,]    5
## [3361,]    9
## [3362,]    4
## [3363,]    5
## [3364,]    7
## [3365,]    4
## [3366,]    6
## [3367,]    6
## [3368,]    7
## [3369,]    6
## [3370,]    6
## [3371,]    7
## [3372,]    9
## [3373,]    9
## [3374,]    8
## [3375,]    8
## [3376,]    9
## [3377,]    4
## [3378,]    6
## [3379,]    5
## [3380,]    6
## [3381,]    8
## [3382,]    5
## [3383,]    4
## [3384,]    7
## [3385,]    5
## [3386,]    3
## [3387,]    5
## [3388,]    4
## [3389,]    5
## [3390,]    3
## [3391,]    8
## [3392,]    7
## [3393,]    3
## [3394,]    8
## [3395,]    8
## [3396,]    8
## [3397,]    7
## [3398,]    9
## [3399,]    8
## [3400,]    8
## [3401,]    5
## [3402,]    6
## [3403,]    8
## [3404,]    8
## [3405,]    7
## [3406,]    9
## [3407,]    5
## [3408,]    5
## [3409,]    7
## [3410,]    5
## [3411,]    7
## [3412,]    5
## [3413,]    6
## [3414,]    4
## [3415,]    2
## [3416,]    3
## [3417,]    3
## [3418,]    9
## [3419,]    4
## [3420,]    7
## [3421,]    8
## [3422,]    6
## [3423,]    7
## [3424,]    7
## [3425,]    7
## [3426,]    8
## [3427,]    7
## [3428,]    9
## [3429,]    7
## [3430,]    9
## [3431,]    9
## [3432,]    7
## [3433,]    9
## [3434,]    7
## [3435,]    8
## [3436,]    8
## [3437,]    4
## [3438,]    5
## [3439,]    5
## [3440,]    4
## [3441,]    6
## [3442,]    5
## [3443,]    4
## [3444,]    7
## [3445,]    6
## [3446,]    7
## [3447,]    7
## [3448,]    1
## [3449,]    5
## [3450,]    5
## [3451,]    6
## [3452,]    7
## [3453,]    7
## [3454,]    7
## [3455,]    8
## [3456,]    8
## [3457,]    8
## [3458,]    7
## [3459,]    7
## [3460,]    5
## [3461,]    9
## [3462,]    7
## [3463,]    9
## [3464,]    1
## [3465,]    2
## [3466,]    3
## [3467,]    4
## [3468,]    2
## [3469,]    5
## [3470,]    3
## [3471,]    3
## [3472,]    3
## [3473,]    5
## [3474,]    5
## [3475,]    7
## [3476,]    8
## [3477,]    8
## [3478,]    8
## [3479,]    3
## [3480,]    5
## [3481,]    5
## [3482,]    6
## [3483,]    8
## [3484,]    5
## [3485,]    3
## [3486,]    8
## [3487,]    6
## [3488,]    7
## [3489,]    7
## [3490,]    6
## [3491,]    9
## [3492,]    2
## [3493,]    5
## [3494,]    4
## [3495,]    6
## [3496,]    5
## [3497,]    5
## [3498,]    3
## [3499,]    3
## [3500,]    7
## [3501,]    7
## [3502,]    7
## [3503,]    7
## [3504,]    6
## [3505,]    9
## [3506,]    6
## [3507,]    6
## [3508,]    4
## [3509,]    8
## [3510,]    5
## [3511,]    9
## [3512,]    5
## [3513,]    8
## [3514,]    7
## [3515,]    4
## [3516,]    7
## [3517,]    5
## [3518,]    4
## [3519,]    4
## [3520,]    4
## [3521,]    7
## [3522,]    5
## [3523,]    7
## [3524,]    8
## [3525,]    8
## [3526,]    7
## [3527,]    2
## [3528,]    5
## [3529,]    4
## [3530,]    6
## [3531,]    5
## [3532,]    2
## [3533,]    4
## [3534,]    7
## [3535,]    9
## [3536,]    6
## [3537,]    8
## [3538,]    6
## [3539,]    6
## [3540,]    8
## [3541,]    5
## [3542,]    8
## [3543,]    3
## [3544,]    5
## [3545,]    6
## [3546,]    2
## [3547,]    5
## [3548,]    8
## [3549,]    7
## [3550,]    8
## [3551,]    9
## [3552,]    8
## [3553,]    7
## [3554,]    7
## [3555,]    7
## [3556,]    5
## [3557,]    7
## [3558,]    6
## [3559,]    8
## [3560,]    6
## [3561,]    8
## [3562,]    8
## [3563,]    9
## [3564,]    7
## [3565,]    8
## [3566,]    5
## [3567,]    5
## [3568,]    5
## [3569,]    4
## [3570,]    7
## [3571,]    6
## [3572,]    7
## [3573,]    4
## [3574,]    8
## [3575,]    8
## [3576,]    7
## [3577,]    7
## [3578,]    3
## [3579,]    7
## [3580,]    8
## [3581,]    8
## [3582,]    5
## [3583,]    8
## [3584,]    8
## [3585,]    4
## [3586,]    8
## [3587,]    5
## [3588,]    6
## [3589,]    6
## [3590,]    3
## [3591,]    5
## [3592,]    4
## [3593,]    7
## [3594,]    5
## [3595,]    2
## [3596,]    3
## [3597,]    7
## [3598,]    8
## [3599,]    7
## [3600,]    3
## [3601,]    9
## [3602,]    4
## [3603,]    6
## [3604,]    6
## [3605,]    8
## [3606,]    8
## [3607,]    2
## [3608,]    4
## [3609,]    6
## [3610,]    6
## [3611,]    2
## [3612,]    8
## [3613,]    4
## [3614,]    5
## [3615,]    5
## [3616,]    1
## [3617,]    7
## [3618,]    5
## [3619,]    4
## [3620,]    6
## [3621,]    7
## [3622,]    7
## [3623,]    5
## [3624,]    4
## [3625,]    8
## [3626,]    5
## [3627,]    7
## [3628,]    9
## [3629,]    8
## [3630,]    8
## [3631,]    3
## [3632,]    6
## [3633,]    7
## [3634,]    6
## [3635,]    4
## [3636,]    5
## [3637,]    5
## [3638,]    6
## [3639,]    4
## [3640,]    6
## [3641,]    4
## [3642,]    7
## [3643,]    5
## [3644,]    5
## [3645,]    6
## [3646,]    6
## [3647,]    6
## [3648,]    7
## [3649,]    7
## [3650,]    7
## [3651,]    7
## [3652,]    9
## [3653,]    7
## [3654,]    8
## [3655,]    8
## [3656,]    8
## [3657,]    8
## [3658,]    8
## [3659,]    5
## [3660,]    5
## [3661,]    6
## [3662,]    3
## [3663,]    8
## [3664,]    3
## [3665,]    4
## [3666,]    6
## [3667,]    5
## [3668,]    7
## [3669,]    7
## [3670,]    8
## [3671,]    5
## [3672,]    4
## [3673,]    4
## [3674,]    5
## [3675,]    8
## [3676,]    3
## [3677,]    6
## [3678,]    9
## [3679,]    5
## [3680,]    1
## [3681,]    9
## [3682,]    7
## [3683,]    6
## [3684,]    4
## [3685,]    7
## [3686,]    6
## [3687,]    8
## [3688,]    5
## [3689,]    6
## [3690,]    5
## [3691,]    9
## [3692,]    5
## [3693,]    6
## [3694,]    7
## [3695,]    7
## [3696,]    8
## [3697,]    6
## [3698,]    4
## [3699,]    6
## [3700,]    7
## [3701,]    6
## [3702,]    4
## [3703,]    5
## [3704,]    5
## [3705,]    5
## [3706,]    6
## [3707,]    7
## [3708,]    3
## [3709,]    6
## [3710,]    8
## [3711,]    9
## [3712,]    9
## [3713,]    0
## [3714,]    2
## [3715,]    8
## [3716,]    6
## [3717,]    7
## [3718,]    9
## [3719,]    6
## [3720,]    7
## [3721,]    9
## [3722,]    5
## [3723,]    7
## [3724,]    5
## [3725,]    8
## [3726,]    4
## [3727,]    5
## [3728,]    4
## [3729,]    8
## [3730,]    9
## [3731,]    7
## [3732,]    6
## [3733,]    8
## [3734,]    9
## [3735,]    4
## [3736,]    5
## [3737,]    6
## [3738,]    7
## [3739,]    6
## [3740,]    4
## [3741,]    5
## [3742,]    6
## [3743,]    6
## [3744,]    3
## [3745,]    2
## [3746,]    6
## [3747,]    5
## [3748,]    7
## [3749,]    6
## [3750,]    9
## [3751,]    9
## [3752,]    7
## [3753,]    6
## [3754,]    7
## [3755,]    8
## [3756,]    7
## [3757,]    9
## [3758,]    7
## [3759,]    7
## [3760,]    7
## [3761,]    6
## [3762,]    6
## [3763,]    8
## [3764,]    5
## [3765,]    7
## [3766,]    5
## [3767,]    7
## [3768,]    6
## [3769,]    4
## [3770,]    4
## [3771,]    8
## [3772,]    7
## [3773,]    4
## [3774,]    5
## [3775,]    7
## [3776,]    9
## [3777,]    6
## [3778,]    6
## [3779,]    9
## [3780,]    9
## [3781,]    6
## [3782,]    7
## [3783,]    3
## [3784,]    4
## [3785,]    2
## [3786,]    2
## [3787,]    4
## [3788,]    4
## [3789,]    4
## [3790,]    5
## [3791,]    3
## [3792,]    5
## [3793,]    6
## [3794,]    6
## [3795,]    6
## [3796,]    5
## [3797,]    5
## [3798,]    6
## [3799,]    5
## [3800,]    4
## [3801,]    8
## [3802,]    8
## [3803,]    8
## [3804,]    7
## [3805,]    6
## [3806,]    3
## [3807,]    6
## [3808,]    2
## [3809,]    7
## [3810,]    3
## [3811,]    3
## [3812,]    4
## [3813,]    6
## [3814,]    7
## [3815,]    6
## [3816,]    3
## [3817,]    6
## [3818,]    5
## [3819,]    5
## [3820,]    3
## [3821,]    5
## [3822,]    8
## [3823,]    0
## [3824,]    7
## [3825,]    7
## [3826,]    8
## [3827,]    4
## [3828,]    6
## [3829,]    5
## [3830,]    3
## [3831,]    3
## [3832,]    3
## [3833,]    2
## [3834,]    7
## [3835,]    7
## [3836,]    7
## [3837,]    5
## [3838,]    8
## [3839,]    6
## [3840,]    5
## [3841,]    3
## [3842,]    7
## [3843,]    6
## [3844,]    4
## [3845,]    6
## [3846,]    8
## [3847,]    6
## [3848,]    7
## [3849,]    7
## [3850,]    5
## [3851,]    5
## [3852,]    2
## [3853,]    5
## [3854,]    5
## [3855,]    6
## [3856,]    8
## [3857,]    8
## [3858,]    8
## [3859,]    9
## [3860,]    8
## [3861,]    7
## [3862,]    9
## [3863,]    4
## [3864,]    7
## [3865,]    7
## [3866,]    7
## [3867,]    6
## [3868,]    8
## [3869,]    8
## [3870,]    6
## [3871,]    7
## [3872,]    8
## [3873,]    6
## [3874,]    8
## [3875,]    5
## [3876,]    5
## [3877,]    7
## [3878,]    9
## [3879,]    6
## [3880,]    5
## [3881,]    6
## [3882,]    6
## [3883,]    4
## [3884,]    2
## [3885,]    6
## [3886,]    9
## [3887,]    6
## [3888,]    8
## [3889,]    8
## [3890,]    6
## [3891,]    7
## [3892,]    8
## [3893,]    6
## [3894,]    6
## [3895,]    7
## [3896,]    5
## [3897,]    3
## [3898,]    7
## [3899,]    7
## [3900,]    7
## [3901,]    9
## [3902,]    6
## [3903,]    6
## [3904,]    7
## [3905,]    7
## [3906,]    6
## [3907,]    2
## [3908,]    7
## [3909,]    5
## [3910,]    3
## [3911,]    7
## [3912,]    6
## [3913,]    7
## [3914,]    6
## [3915,]    7
## [3916,]    7
## [3917,]    6
## [3918,]    8
## [3919,]    7
## [3920,]    5
## [3921,]    7
## [3922,]    7
## [3923,]    5
## [3924,]    7
## [3925,]    8
## [3926,]    6
## [3927,]    4
## [3928,]    6
## [3929,]    7
## [3930,]    4
## [3931,]    4
## [3932,]    3
## [3933,]    7
## [3934,]    6
## [3935,]    8
## [3936,]    5
## [3937,]    4
## [3938,]    6
## [3939,]    7
## [3940,]    5
## [3941,]    6
## [3942,]    3
## [3943,]    5
## [3944,]    3
## [3945,]    5
## [3946,]    7
## [3947,]    4
## [3948,]    8
## [3949,]    2
## [3950,]    4
## [3951,]    7
## [3952,]    8
## [3953,]    7
## [3954,]    5
## [3955,]    3
## [3956,]    6
## [3957,]    7
## [3958,]    4
## [3959,]    3
## [3960,]    8
## [3961,]    5
## [3962,]    7
## [3963,]    7
## [3964,]    8
## [3965,]    6
## [3966,]    5
## [3967,]    3
## [3968,]    9
## [3969,]    5
## [3970,]    2
## [3971,]    2
## [3972,]    4
## [3973,]    1
## [3974,]    4
## [3975,]    5
## [3976,]    8
## [3977,]    4
## [3978,]    4
## [3979,]    5
## [3980,]    9
## [3981,]    4
## [3982,]    5
## [3983,]    5
## [3984,]    7
## [3985,]    4
## [3986,]    9
## [3987,]    5
## [3988,]    7
## [3989,]    7
## [3990,]    4
## [3991,]    7
## [3992,]    0
## [3993,]    4
## [3994,]    5
## [3995,]    5
## [3996,]    7
## [3997,]    6
## [3998,]    9
## [3999,]    8
## [4000,]    7

data.frame(obs = ppd) %>% 
  ggplot(aes(x=obs)) +
  geom_histogram() +
  labs(x="Observed water (out of 9)")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.


Sampling parameters from the prior

Oops, we’ve jumped ahead of ourselves! Best practice is to simulate values from your prior first and check to see if those priors are reasonable.

m1p <-
  brm(data = list(w = 6),                            
      family = binomial(link = "identity"),          
      w | trials(9) ~ 0 + Intercept,                 
      prior(beta(1, 1), class = b, lb = 0, ub = 1),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,          
      sample_prior = "only")
samples_from_prior = as_draws_df(m1p)
samples_from_prior %>% 
  ggplot(aes(x=b_Intercept)) +
  geom_density(fill = "grey", color = "white") +
  labs(x="Proportion water", title="Prior")


Simulating observations from the prior

We may also want the PRIOR PREDICTIVE DISTRIBUTION which is the expected observiations given our prior.

prior_pd = posterior_predict(m1p)
data.frame(obs = prior_pd) %>% 
  ggplot(aes(x=obs)) +
  geom_histogram() +
  labs(x="Observed water (out of 9)")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Simulating from your priors – prior predictive simulation – is an essential part of modeling. This allows you to see what your choices imply about the data. You’ll be able to diagnose bad choices.


an aside about learning in R

At this point in the course, I’m going to start throwing a lot of code at you. Do I expect you to memorize this code? Of course not.

Do you need to understand every single thing that’s happening in the code? Nope.

But, you’ll learn a lot by taking the time to figure out what’s happening in a code chunk. Class time will frequently include exercises where I ask you to adapt code I’ve shared in the slides to a new dataset or to answer a new problem. When doing so, go back through the old code and figure out what’s going on. Run the code one line at a time. Always observe the output and take some time to look at the object that was created or modified. Here are some functions that will be extremely useful:

str() # what kind of object is this? what is its structure?
dim() # what are the dimensions (rows/columns) of this object
head() # give me the first bit of this object
str(prior_pd)
##  int [1:4000, 1] 6 0 1 1 9 7 0 5 4 5 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : NULL
##   ..$ : NULL
dim(prior_pd)
## [1] 4000    1
head(prior_pd)
##      [,1]
## [1,]    6
## [2,]    0
## [3,]    1
## [4,]    1
## [5,]    9
## [6,]    7

Continous outcomes

The globe tossing example is cute and easy to work with, but let’s move towards the kinds of variables we more frequently work with. Let’s create a model for some outcome, \(y\) that is continuous.

\[\begin{align*} y_i &\sim \text{Normal}(\mu, \sigma) \\ \mu &\sim \text{Normal}(0, 10) \\ \sigma &\sim \text{Uniform}(0, 5) \end{align*}\]

set.seed(9)
y = rnorm(n = 31, mean = 4, sd = .5)
m2 = brm(
  data = list(y=y),
  family = gaussian,
  y ~ 1,
  prior = c(prior( normal(0,10), class=Intercept),
            prior( uniform(0,5), class=sigma)),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,
  file = here("files/models/m21.2")
)

An example: weight and height

Using the Howell data (don’t load the rethinking package because it can interfere with brms).

data("Howell1", package = "rethinking")
d <- Howell1
str(d)
## 'data.frame':    544 obs. of  4 variables:
##  $ height: num  152 140 137 157 145 ...
##  $ weight: num  47.8 36.5 31.9 53 41.3 ...
##  $ age   : num  63 63 65 41 51 35 32 27 19 54 ...
##  $ male  : int  1 0 0 1 0 1 0 1 0 1 ...
library(measurements)
d$height <- conv_unit(d$height, from = "cm", to = "feet")
d$weight <- conv_unit(d$weight, from = "kg", to = "lbs")
rethinking::precis(d)
##              mean         sd      5.5%    94.5%      histogram
## height  4.5362072  0.9055921  2.661042   5.4375      ▁▁▁▁▂▂▇▇▁
## weight 78.5079631 32.4502290 20.636856 120.1583 ▁▂▃▂▂▁▁▃▅▇▇▃▂▁
## age    29.3443934 20.7468882  1.000000  66.1350      ▇▅▅▃▅▂▂▁▁
## male    0.4724265  0.4996986  0.000000   1.0000     ▇▁▁▁▁▁▁▁▁▇
d2 <- d[ d$age >= 18, ]

exercise

Write a mathematical model for the weights in this data set. (Don’t worry about predicting from other variables yet.)

solution

\[\begin{align*} w &\sim \text{Normal}(\mu, \sigma) \\ \mu &\sim \text{Normal}(130, 20) \\ \sigma &\sim \text{Uniform}(0, 25) \\ \end{align*}\]


\[\begin{align*} w &\sim \text{Normal}(\mu, \sigma) \\ \mu &\sim \text{Normal}(130, 20) \\ \sigma &\sim \text{Uniform}(0, 25) \\ \end{align*}\]

exercise

Simulate from your priors (parameters values and prior predictive values).

solution

Sample from your priors:

m3p = brm(
  data = d2,
  family = gaussian,
  weight ~ 1,
  prior = c(prior( normal(130,20), class=Intercept),
            prior( uniform(0,25), class=sigma, lb=0, ub=25)),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,
  sample_prior = "only")

Sampling parameter estimates for the Intercept.

pairs(m3p)
## Warning: Only one chain in 'x'. This plot is more useful with multiple chains.

This is a different (shorter) way to plot your posterior. Good things: it automatically includes the scatterplot so you can see the implications of how these parameters correlate. Bad things: not customizable and not useable when you have a lot of parameters.


Simulate values of weight.

prior_pd = posterior_predict(m3p)
dim(prior_pd)
## [1] 4000  352
as.data.frame(prior_pd) %>% 
  pivot_longer(everything()) %>% 
  ggplot(aes(x=value)) +
  geom_histogram() +
  labs(x="Expected observed weights (based on prior)")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.


Another shorter way:

pp_check(m3p)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.


Fit the model

m3 = brm(
  data = d2,
  family = gaussian,
  weight ~ 1,
  prior = c(prior( normal(130,20), class=Intercept),
            prior( uniform(0,25), class=sigma, lb=0, ub=25)),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,
  file = here("files/models/m21.3"))
posterior_summary(m3)
##                 Estimate  Est.Error         Q2.5        Q97.5
## b_Intercept    99.221909 0.75720356    97.751074   100.737488
## sigma          14.291987 0.55198322    13.254363    15.428183
## Intercept      99.221909 0.75720356    97.751074   100.737488
## lprior         -8.318377 0.05827127    -8.433538    -8.203915
## lp__        -1441.294276 1.02695551 -1444.235605 -1440.292326

pairs(m3)
## Warning: Only one chain in 'x'. This plot is more useful with multiple chains.


posterior_predict(m3) %>% 
  as.data.frame() %>% 
  pivot_longer(everything()) %>% 
  ggplot(aes(x=value)) +
  geom_density(fill = "grey", color = "white") +
  geom_density( aes(x = weight), data=d2, inherit.aes = F) 


pp_check(m3)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.


Adding in a linear component

We might assume that height and weight are associated with each other. Indeed, within our sample:

plot(d2$weight ~ d2$height)


exercise

Update your mathematical model to incorporate height.

\[\begin{align*} w_i &\sim \text{Normal}(\mu_i, \sigma) \\ \mu_i &= \alpha + \beta h_i \\ \alpha &\sim \text{Normal}(??, ??) \\ \beta &\sim \text{Normal}(0, 25) \\ \sigma &\sim \text{Uniform}(0, 25) \\ \end{align*}\]

\(=\) is deterministic – once we know other variables, \(\mu_i\) is known with certainty

made-up parameters are the targets of learning


exercise

Update your mathematical model to incorporate height.

\[\begin{align*} w_i &\sim \text{Normal}(\mu_i, \sigma) \\ \mu_i &= \alpha + \beta (h_i - \bar{h}) \\ \alpha &\sim \text{Normal}(130, 20) \\ \beta &\sim \text{Normal}(0, 25) \\ \sigma &\sim \text{Uniform}(0, 25) \\ \end{align*}\]

\(=\) is deterministic – once we know other variables, \(\mu_i\) is known with certainty

made-up parameters are the targets of learning


\[\begin{align*} w_i &\sim \text{Normal}(\mu_i, \sigma) \\ \mu_i &= \alpha + \beta (h_i - \bar{h}) \\ \alpha &\sim \text{Normal}(130, 20) \\ \beta &\sim \text{Normal}(0, 25) \\ \sigma &\sim \text{Uniform}(0, 25) \\ \end{align*}\]

To update our brms code:

d2$height_c = d2$height - mean(d2$height)
  
m4p = brm(
  data = d2,
  family = gaussian,
  weight ~ 1 + height_c,
  prior = c(prior( normal(130,20), class=Intercept),
            prior( normal(0,25),   class=b),
            prior( uniform(0,25),  class=sigma, lb=0, ub=25)),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,
  sample_prior = "only")

set.seed(9)
samples_from_prior = as_draws_df(m4p)
str(samples_from_prior)
## draws_df [4,000 × 9] (S3: draws_df/draws/tbl_df/tbl/data.frame)
##  $ b_Intercept: num [1:4000] 167 155.6 149.6 88.7 176.4 ...
##  $ b_height_c : num [1:4000] -57.9 -57.2 -18.5 52 -47.3 ...
##  $ sigma      : num [1:4000] 14.1 11.8 13.3 2.1 21.1 ...
##  $ Intercept  : num [1:4000] 167 155.6 149.6 88.7 176.4 ...
##  $ lprior     : num [1:4000] -15.7 -14.7 -12 -15.6 -15.8 ...
##  $ lp__       : num [1:4000] -13.8 -12.9 -10.2 -14.9 -14.6 ...
##  $ .chain     : int [1:4000] 1 1 1 1 1 1 1 1 1 1 ...
##  $ .iteration : int [1:4000] 1 2 3 4 5 6 7 8 9 10 ...
##  $ .draw      : int [1:4000] 1 2 3 4 5 6 7 8 9 10 ...

d2 %>% 
  ggplot(aes(x=height_c, y=weight)) +
  geom_blank() +
  geom_abline( aes(intercept=b_Intercept, slope=b_height_c), 
               data=samples_from_prior[1:50, ],
               alpha=.3) +
  scale_x_continuous(name = "height(feet)", 
                     breaks=seq(4,6,by=.5)-mean(d2$height),
                     labels=seq(4,6,by=.5))

Describe in words what’s wrong with our priors.

Slope should not be negative. How can we fix this?

Could use a uniform distribution bounded by 0.


m4p = brm(
  data = d2,
  family = gaussian,
  weight ~ 1 + height_c,
  prior = c(prior( normal(130,20), class=Intercept),
            prior( uniform(0,25),   class=b),
            prior( uniform(0,25),  class=sigma, lb=0, ub=25)),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,
  sample_prior = "only")


exercise

Fit the new weight model to the data.

solution

m4 = brm(
  data = d2,
  family = gaussian,
  weight ~ 1 + height_c,
  prior = c(prior( normal(130,20), class=Intercept),
            prior( uniform(0,25),   class=b),
            prior( uniform(0,25),  class=sigma, lb=0, ub=25)),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,
  file=here("files/models/m21.4"))
posterior_summary(m4)
##                Estimate Est.Error        Q2.5       Q97.5
## b_Intercept    99.21508 0.5154139    98.17317   100.20955
## b_height_c     24.74661 0.2609169    24.05150    24.99355
## sigma          10.36090 0.3898531     9.65234    11.12228
## Intercept      99.21508 0.5154139    98.17317   100.20955
## lprior        -11.53739 0.0396881   -11.61861   -11.46176
## lp__        -1332.15882 1.3936194 -1335.90639 -1330.52002

exercise

Draw lines from the posterior distribution and plot with the data.

solution

set.seed(9)
samples_from_posterior = as_draws_df(m4)
d2%>% 
  ggplot(aes(x=height_c, y=weight)) +
  geom_point(size=.5) +
  geom_abline( aes(intercept=b_Intercept, slope=b_height_c), 
               data=samples_from_posterior[1:50, ],
               alpha=.3,
               color="#1c5253") +
  scale_x_continuous(name = "height(feet)", 
                     breaks=seq(4,6,by=.5)-mean(d2$height),
                     labels=seq(4,6,by=.5))


A side note: a major concern or critique of Bayesian analysis is that the subjectivity of the priors allow for nefarious behavior. “Putting our thumbs on the scale,” so to speak. But priors are quickly overwhelmed by data. Case in point:

m4e = brm(
  data = d2,
  family = gaussian,
  weight ~ 1 + height_c,
  prior = c(prior( normal(130,20), class=Intercept),
            prior( normal(-5,5),   class=b),
            prior( uniform(0,25),  class=sigma, lb=0, ub=25)),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,
  file=here("files/models/m21.4e"))
posterior_summary(m4e)
##                 Estimate Est.Error         Q2.5       Q97.5
## b_Intercept    99.197733 0.5035653    98.214092   100.15932
## b_height_c     35.775818 1.8832462    32.177572    39.50070
## sigma           9.529429 0.3687178     8.839379    10.27233
## Intercept      99.197733 0.5035653    98.214092   100.15932
## lprior        -44.172476 3.0738968   -50.456409   -38.49994
## lp__        -1334.629214 1.1849503 -1337.645046 -1333.23509

You’ll only really get into trouble with uniform priors that have a boundary, if true population parameter is outside your boundary. A good rule of thumb is to avoid the uniform distribution. We’ll cover other options for priors for \(\sigma\) in future lectures, but as a preview, the exponential distribution works very well for this!